text:
module examples_fortran
  use, intrinsic :: iso_c_binding
  use, intrinsic :: omp_lib, only : omp_get_num_threads, omp_get_thread_num
  implicit none

  contains

  subroutine f90_hello() &
  bind(C, name="f90_hello")
    integer :: tid, nthreads

    !$omp parallel private(tid)
    nthreads = omp_get_num_threads()
    tid = omp_get_thread_num()

    call intpr("Hello from thread", -1, tid, 1)
    !$omp end parallel
  end subroutine

  function f90_sum(x, xlen) result(s) &
  bind(C, name="f90_sum")
    real(kind=c_double) :: s
    integer(kind=c_int), intent(in), value :: xlen
    real(kind=c_double), intent(in) :: x(xlen)
    integer i

    s = 0.0d0

    !$omp parallel do default(shared) private(i) reduction(+:s)
    do i=1, xlen
      s = s + x(i)
    end do
  end function

  subroutine f90_sweep(m, n, x, vec, ret) &
  bind(C, name="f90_sweep")
    integer(kind=c_int), intent(in), value :: m, n
    real(kind=c_double), intent(in) :: x(m, n), vec(m)
    real(kind=c_double), intent(out) :: ret(m, n)
    integer :: i, j

    !$omp parallel do default(shared) private(i, j)
    do j=1, n
      do i=1, m
        ret(i, j) = x(i, j) - vec(i)
      end do
    end do
    !$omp end parallel do
  end subroutine

  function f90_primesbelow(n) result(nprimes) &
  bind(C, name="f90_primesbelow")
    integer(kind=c_int) :: nprimes
    integer(kind=c_int), intent(in), value :: n
    integer i, j, isprime

    nprimes = 1

    !$omp parallel do private(i, j, isprime) reduction(+:nprimes)
    do i=3, n, 2
      isprime = 1

      do j=3, i-1, 2
        if (mod(i, j) == 0) then
          isprime = 0
          exit
        end if
      end do

      if (isprime == 1) nprimes = nprimes+1
    end do
  end function

end module

max_stars_repo_name: wrathematics/Romp

max_issues_repo_name: wrathematics/Romp

max_forks_repo_name: wrathematics/Romp


!************************************************************************
!                  Fortran 95 OPA Nesting tools			    			*
!					                                    				*
!     Copyright (C) 2005 Florian Lemari�(Florian.Lemarie@imag.fr)	    *
!									                                    *
!************************************************************************
!
!********************************************************************************
!		                                        								*
!                              module io_netcdf		    						*
!						                                        				*
!                  NetCDF Fortran 90 read/write interface	    				*
!                  using input/output functions provided						*
!                                 by unidata									*
!						                                        				*
!http://my.unidata.ucar.edu/content/software/netcdf/docs/netcdf-f90/index.html	*
!									                                        	*
!********************************************************************************
!
!
!
MODULE io_netcdf
  !
  USE netcdf
  USE types
  !
  INTERFACE Read_Ncdf_var
     MODULE PROCEDURE Read_Ncdf_var1d_Real,       &
                      Read_Ncdf_var2d_Real,       &
                      Read_Ncdf_var2d_Real_bis,   &
                      Read_Ncdf_var3d_Real,       &
                      Read_Ncdf_var4d_Real,       &
                      Read_Ncdf_var3d_Real_t,     &
                      Read_Ncdf_var4d_Real_t,     &
                      Read_Ncdf_var4d_Real_nt,    &
                      Read_Ncdf_var1d_Int,        &
                      Read_Ncdf_var2d_Int,        &
                      Read_Ncdf_var3d_Int,        &
                      Read_Ncdf_var4d_Int,        &
                      Read_Ncdf_var0d_Int,        &
                      Read_Ncdf_var0d_Real
  END INTERFACE
  !
  INTERFACE Write_Ncdf_var
     MODULE PROCEDURE Write_Ncdf_var1d_Real,      &
					  Write_Ncdf_var2d_Real,      &
					  Write_Ncdf_var3d_Real,      &
					  Write_Ncdf_var4d_Real,      &
					  Write_Ncdf_var3d_Real_t,    &
					  Write_Ncdf_var4d_Real_t,    &
					  Write_Ncdf_var4d_Real_nt,   &
					  Write_Ncdf_var2d_Real_bis,  &
					  Write_Ncdf_var1d_Int,       &
					  Write_Ncdf_var2d_Int,       &
					  Write_Ncdf_var3d_Int,       &
					  Write_Ncdf_var4d_Int,       &
					  Write_Ncdf_var0d_Real
  END INTERFACE
  !
  INTERFACE Copy_Ncdf_att
     MODULE PROCEDURE Copy_Ncdf_att_latlon,Copy_Ncdf_att_var
  END INTERFACE
  !
CONTAINS
  !
  !****************************************************************
  !                 subroutine Read_Ncdf_dim					*
  !				                                				*
  !     subroutine to retrieve value of a given dimension 		*
  !	                                							*
  !     dimname : name of dimension to retrieve		        	*
  !     file    : netcdf file name				                *
  !     dimval  : value of the required dimension   			*
  !								                                *
  !****************************************************************
  !
  SUBROUTINE Read_Ncdf_dim(dimname,file,dimval)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: dimname,file
    INTEGER    :: dimval
    INTEGER ncid,status,dimid
    !
    status = nf90_open(file,NF90_NOWRITE,ncid)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_dimid(ncid,dimname,dimid)
    status = nf90_inquire_dimension(ncid,dimid,len=dimval)
    status = nf90_close(ncid)
    !
  END SUBROUTINE Read_Ncdf_dim
  !
  !
  !
  !****************************************************************
  !                subroutine Write_Ncdf_dim					*
  !		                                 						*
  !      subroutine to write a dimension in a given file 		*
  !		                                						*
  !     dimname : name of dimension to initialize    			*
  !     file    : netcdf file name				                *
  !     dimval  : value of the dimension to write	     		*
  !								                                *
  !****************************************************************
  SUBROUTINE Write_Ncdf_dim(dimname,file,dimval)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: dimname,file
    INTEGER    :: dimval
    !
    ! local variables
    !
    INTEGER ncid,status,dimid
    !
    status = nf90_open(file,NF90_WRITE,ncid)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_redef(ncid)
    IF(dimval.EQ.0) THEN
       status = nf90_def_dim(ncid,dimname,nf90_unlimited,dimid)
    ELSE
       status = nf90_def_dim(ncid,dimname,dimval,dimid)
    END IF
    !
    status = nf90_enddef(ncid)
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Write_Ncdf_dim
  !
  !
  !
  !****************************************************************
  !                 subroutine Read_Ncdf_var					   *
  !		                                						   *
  !     subroutine to retrieve values of a given variable 		   *
  !				                                				   *
  !     varname : name of variable to retrieve	        		   *
  !     file    : netcdf file name			                	   *
  !     tabvar  : array containing values of the required variable *
  !								                                   *
  !****************************************************************
  SUBROUTINE Read_Ncdf_var1d_Real(varname,file,tabvar)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file
    REAL*8, DIMENSION(:), POINTER :: tabvar
    !
    !local variables
    !
    INTEGER, DIMENSION(1) :: dimID
    INTEGER :: dim1
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_NOWRITE,ncid)
    !
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_varid(ncid,varname,varid)
    status=nf90_inquire_variable(ncid,varid,dimids=dimID)
    status=nf90_inquire_dimension(ncid,dimID(1),len=dim1)
    !
    IF(.NOT. ASSOCIATED(tabvar)) THEN
       ALLOCATE(tabvar(dim1))
    ELSE
       IF( ANY(SHAPE(tabvar)/=(/dim1/)) ) THEN
          DEALLOCATE(tabvar)
          ALLOCATE(tabvar(dim1))
       ENDIF
    ENDIF
    !
    status=nf90_get_var(ncid,varid,tabvar)
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Read_Ncdf_var1d_Real
  !
  !
  !
  !**************************************************************
  !         subroutine Read_Ncdf_var2d_real
  !**************************************************************
  SUBROUTINE Read_Ncdf_var2d_Real(varname,file,tabvar)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file
    REAL*8, DIMENSION(:,:), POINTER :: tabvar
    !local variables
    INTEGER, DIMENSION(10) :: dimIDS
    INTEGER :: dim1,dim2
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_NOWRITE,ncid)
    !
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_varid(ncid,varname,varid)
    status=nf90_inquire_variable(ncid,varid,dimids=dimIDS)
    status=nf90_inquire_dimension(ncid,dimIDS(1),len=dim1)
    status=nf90_inquire_dimension(ncid,dimIDS(2),len=dim2)
    !
    IF(.NOT. ASSOCIATED(tabvar)) THEN
       ALLOCATE(tabvar(dim1,dim2))
    ELSE
       IF( ANY(SHAPE(tabvar)/=(/dim1,dim2/)) ) THEN
          DEALLOCATE(tabvar)
          ALLOCATE(tabvar(dim1,dim2))
       ENDIF
    ENDIF
    !
    status=nf90_get_var(ncid,varid,tabvar)
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Read_Ncdf_var2d_Real
  !
  !
  !
  !**************************************************************
  !   subroutine Read_Ncdf_var2d_real_bis
  !**************************************************************
  SUBROUTINE Read_Ncdf_var2d_Real_bis(varname,file,tabvar,strt,cnt)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file
    REAL*8, DIMENSION(:,:), POINTER :: tabvar
    !local variables
    INTEGER, DIMENSION(10) :: dimIDS
    INTEGER, DIMENSION(2) :: strt,cnt
    INTEGER :: dim1,dim2
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_NOWRITE,ncid)
    !
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    dim1 = cnt(1)
    dim2 = cnt(2)
    !
    status = nf90_inq_varid(ncid,varname,varid)
    status=nf90_inquire_variable(ncid,varid,dimids=dimIDS)
    !
    IF(.NOT. ASSOCIATED(tabvar)) THEN
       ALLOCATE(tabvar(dim1,dim2))
    ELSE
       IF( ANY(SHAPE(tabvar)/=(/dim1,dim2/)) ) THEN
          DEALLOCATE(tabvar)
          ALLOCATE(tabvar(dim1,dim2))
       ENDIF
    ENDIF
    !
    status=nf90_get_var(ncid,varid,tabvar,start = strt,count = cnt)
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Read_Ncdf_var2d_Real_bis
  !
  !
  !
  !**************************************************************
  !   subroutine Read_Ncdf_var3d_real
  !**************************************************************
  SUBROUTINE Read_Ncdf_var3d_Real(varname,file,tabvar)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file
    REAL*8, DIMENSION(:,:,:), POINTER :: tabvar
    !
    !local variables
    !
    INTEGER, DIMENSION(10) :: dimIDS
    INTEGER :: dim1,dim2,dim3
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_NOWRITE,ncid)
    !
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_varid(ncid,varname,varid)
    status=nf90_inquire_variable(ncid,varid,dimids=dimIDS)
    status=nf90_inquire_dimension(ncid,dimIDS(1),len=dim1)
    status=nf90_inquire_dimension(ncid,dimIDS(2),len=dim2)
    status=nf90_inquire_dimension(ncid,dimIDS(3),len=dim3)
    !
    IF(.NOT. ASSOCIATED(tabvar)) THEN
       ALLOCATE(tabvar(dim1,dim2,dim3))
    ELSE
       IF( ANY(SHAPE(tabvar) /= (/dim1,dim2,dim3/)) ) THEN
          DEALLOCATE(tabvar)
          ALLOCATE(tabvar(dim1,dim2,dim3))
       ENDIF
    ENDIF
    !
    status=nf90_get_var(ncid,varid,tabvar)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to retrieve netcdf variable : ",TRIM(varname)
       STOP
    ENDIF
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Read_Ncdf_var3d_Real
  !
  !
  !
  !**************************************************************
  !         subroutine Read_Ncdf_var4d_real
  !**************************************************************
  SUBROUTINE Read_Ncdf_var4d_Real(varname,file,tabvar)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file
    REAL*8, DIMENSION(:,:,:,:), POINTER :: tabvar
    !
    !local variables
    !
    INTEGER, DIMENSION(10) :: dimIDS
    INTEGER :: dim1,dim2,dim3,dim4
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_NOWRITE,ncid)
    !
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_varid(ncid,varname,varid)
    status=nf90_inquire_variable(ncid,varid,dimids=dimIDS)
    status=nf90_inquire_dimension(ncid,dimIDS(1),len=dim1)
    status=nf90_inquire_dimension(ncid,dimIDS(2),len=dim2)
    status=nf90_inquire_dimension(ncid,dimIDS(3),len=dim3)
    status=nf90_inquire_dimension(ncid,dimIDS(4),len=dim4)
    !
    IF(.NOT. ASSOCIATED(tabvar)) THEN
       ALLOCATE(tabvar(dim1,dim2,dim3,dim4))
    ELSE
       IF( ANY(SHAPE(tabvar) /= (/dim1,dim2,dim3,dim4/)) ) THEN
          DEALLOCATE(tabvar)
          ALLOCATE(tabvar(dim1,dim2,dim3,dim4))
       ENDIF
    ENDIF
    !
    status=nf90_get_var(ncid,varid,tabvar)
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Read_Ncdf_var4d_Real

  SUBROUTINE Read_Ncdf_var0d_Real(varname,file,tabvar)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file
    REAL*8 :: tabvar
    !
    !local variables
    !
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_NOWRITE,ncid)
    !
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_varid(ncid,varname,varid)
    !
    !
    status=nf90_get_var(ncid,varid,tabvar)
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Read_Ncdf_var0d_Real

  SUBROUTINE Read_Ncdf_var0d_Int(varname,file,tabvar)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file
    INTEGER :: tabvar
    !
    !local variables
    !
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_NOWRITE,ncid)
    !
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_varid(ncid,varname,varid)
    !
    !
    status=nf90_get_var(ncid,varid,tabvar)
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Read_Ncdf_var0d_Int
  !
  !
  !
  !**************************************************************
  !   subroutine Read_Ncdf_var1d_int
  !**************************************************************
  SUBROUTINE Read_Ncdf_var1d_Int(varname,file,tabvar)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file
    INTEGER, DIMENSION(:), POINTER :: tabvar
    !
    !local variables
    !
    INTEGER,DIMENSION(10) :: dimID
    INTEGER :: dim1
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_NOWRITE,ncid)
    !
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_varid(ncid,varname,varid)
    status=nf90_inquire_variable(ncid,varid,dimids=dimID)
    status=nf90_inquire_dimension(ncid,dimID(1),len=dim1)
    !
    IF(.NOT. ASSOCIATED(tabvar)) THEN
       ALLOCATE(tabvar(dim1))
    ELSE
       IF( ANY(SHAPE(tabvar) /= (/dim1/)) ) THEN
          DEALLOCATE(tabvar)
          ALLOCATE(tabvar(dim1))
       ENDIF
    ENDIF
    !
    status=nf90_get_var(ncid,varid,tabvar)
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Read_Ncdf_var1d_Int
  !
  !
  !
  !**************************************************************
  !   subroutine Read_Ncdf_var2d_int
  !**************************************************************
  SUBROUTINE Read_Ncdf_var2d_Int(varname,file,tabvar)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file
    INTEGER, DIMENSION(:,:), POINTER :: tabvar
    !local variables
    INTEGER, DIMENSION(10) :: dimIDS
    INTEGER :: dim1,dim2
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_NOWRITE,ncid)
    !
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_varid(ncid,varname,varid)
    status=nf90_inquire_variable(ncid,varid,dimids=dimIDS)
    status=nf90_inquire_dimension(ncid,dimIDS(1),len=dim1)
    status=nf90_inquire_dimension(ncid,dimIDS(2),len=dim2)
    !
    IF(.NOT. ASSOCIATED(tabvar)) THEN
       ALLOCATE(tabvar(dim1,dim2))
    ELSE
       IF( ANY(SHAPE(tabvar) /= (/dim1,dim2/)) ) THEN
          DEALLOCATE(tabvar)
          ALLOCATE(tabvar(dim1,dim2))
       ENDIF
    ENDIF
    !
    status=nf90_get_var(ncid,varid,tabvar)
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Read_Ncdf_var2d_Int
  !
  !
  !
  !**************************************************************
  !   subroutine Read_Ncdf_var3d_Int
  !**************************************************************
  SUBROUTINE Read_Ncdf_var3d_Int(varname,file,tabvar)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file
    INTEGER, DIMENSION(:,:,:), POINTER :: tabvar
    !
    !local variables
    !
    INTEGER, DIMENSION(10) :: dimIDS
    INTEGER :: dim1,dim2,dim3
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_NOWRITE,ncid)
    !
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_varid(ncid,varname,varid)
    status=nf90_inquire_variable(ncid,varid,dimids=dimIDS)
    status=nf90_inquire_dimension(ncid,dimIDS(1),len=dim1)
    status=nf90_inquire_dimension(ncid,dimIDS(2),len=dim2)
    status=nf90_inquire_dimension(ncid,dimIDS(3),len=dim3)
    !
    IF(.NOT. ASSOCIATED(tabvar)) THEN
       ALLOCATE(tabvar(dim1,dim2,dim3))
    ELSE
       IF( ANY(SHAPE(tabvar) /= (/dim1,dim2,dim3/)) ) THEN
          DEALLOCATE(tabvar)
          ALLOCATE(tabvar(dim1,dim2,dim3))
       ENDIF
    ENDIF
    !
    status=nf90_get_var(ncid,varid,tabvar)
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Read_Ncdf_var3d_Int
  !
  !
  !
  !**************************************************************
  !   subroutine Read_Ncdf_var4d_int
  !**************************************************************
  SUBROUTINE Read_Ncdf_var4d_Int(varname,file,tabvar)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file
    INTEGER, DIMENSION(:,:,:,:), POINTER :: tabvar
    !
    !local variables
    !
    INTEGER, DIMENSION(10) :: dimIDS
    INTEGER :: dim1,dim2,dim3,dim4
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_NOWRITE,ncid)
    !
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_varid(ncid,varname,varid)
    status=nf90_inquire_variable(ncid,varid,dimids=dimIDS)
    status=nf90_inquire_dimension(ncid,dimIDS(1),len=dim1)
    status=nf90_inquire_dimension(ncid,dimIDS(2),len=dim2)
    status=nf90_inquire_dimension(ncid,dimIDS(3),len=dim3)
    status=nf90_inquire_dimension(ncid,dimIDS(4),len=dim4)
    !
    IF(.NOT. ASSOCIATED(tabvar)) THEN
       ALLOCATE(tabvar(dim1,dim2,dim3,dim4))
    ELSE
       IF( ANY(SHAPE(tabvar) /= (/dim1,dim2,dim3,dim4/)) ) THEN
          DEALLOCATE(tabvar)
          ALLOCATE(tabvar(dim1,dim2,dim3,dim4))
       ENDIF
    ENDIF
    !
    status=nf90_get_var(ncid,varid,tabvar)
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Read_Ncdf_var4d_Int
  !
  !
  !
  !****************************************************************
  !                subroutine Write_Ncdf_var					*
  !			                                 					*
  !       subroutine to write a variable in a given file 		*
  !	                                							*
  !     varname : name of variable to store		            	*
  !     dimname : name of dimensions of the given variable	    *
  !     file    : netcdf file name				                *
  !     tabvar  : values of the variable to write		     	*
  !								                                *
  !****************************************************************
  !
  !
  !**************************************************************
  !      subroutine Write_Ncdf_var1d_real
  !**************************************************************
  SUBROUTINE Write_Ncdf_var1d_Real(varname,dimname,file,tabvar,typevar)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file,dimname,typevar
    REAL*8, DIMENSION(:), POINTER :: tabvar
    !
    ! local variables
    !
    INTEGER :: dimid
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_WRITE,ncid)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_dimid(ncid,dimname, dimid)
    status = nf90_inq_varid(ncid,varname,varid)
    status = nf90_redef(ncid)
    SELECT CASE(TRIM(typevar))
    CASE('double')
       status = nf90_def_var(ncid,varname,nf90_double,(/dimid/),varid)
    CASE('float')
       status = nf90_def_var(ncid,varname,nf90_float,(/dimid/),varid)
    END SELECT
    status = nf90_enddef(ncid)
    status = nf90_put_var(ncid,varid,tabvar)
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Write_Ncdf_var1d_Real
  !
  !
  !
  !**************************************************************
  !   subroutine Write_Ncdf_var2d_real
  !**************************************************************
  SUBROUTINE Write_Ncdf_var2d_Real_bis(varname,dimname,file,tabvar,nbdim,typevar)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file,typevar
    INTEGER,INTENT(in) :: nbdim
    CHARACTER(*), DIMENSION(4) :: dimname
    REAL*8, DIMENSION(:,:) :: tabvar
    REAL*8, DIMENSION(:,:,:),POINTER :: tabtemp3d
    REAL*8, DIMENSION(:,:,:,:),POINTER :: tabtemp4d
    !
    ! local variables
    !
    INTEGER :: dimid1,dimid2,dimid3,dimid4
    INTEGER :: status,ncid,ncid2
    INTEGER :: varid,varid2
    !
    IF(nbdim==4) THEN
       ALLOCATE(tabtemp4d(SIZE(tabvar,1),SIZE(tabvar,2),1,1))
       tabtemp4d(:,:,1,1) = tabvar(:,:)
    ELSE IF(nbdim==3) THEN
       ALLOCATE(tabtemp3d(SIZE(tabvar,1),SIZE(tabvar,2),1))
       tabtemp3d(:,:,1) = tabvar(:,:)
    END IF
    !
    status = nf90_open(file,NF90_WRITE,ncid)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_dimid(ncid,dimname(1), dimid1)
    status = nf90_inq_dimid(ncid,dimname(2), dimid2)
    status = nf90_inq_dimid(ncid,dimname(3), dimid3)
    !
    IF(nbdim==4) status = nf90_inq_dimid(ncid,dimname(4), dimid4)
    !
    status = nf90_inq_varid(ncid,varname,varid)
    status = nf90_redef(ncid)
    IF(nbdim==4 .AND. typevar == 'double') THEN
       status = nf90_def_var(ncid,varname,nf90_double,     &
            (/dimid1,dimid2,dimid3,dimid4/),varid)
       !
    ELSE IF(nbdim==4 .AND. typevar == 'float') THEN
       status = nf90_def_var(ncid,varname,nf90_float,     &
            (/dimid1,dimid2,dimid3,dimid4/),varid)
       !
    ELSE IF(nbdim==3 .AND. typevar == 'float') THEN
       status = nf90_def_var(ncid,varname,nf90_float,     &
            (/dimid1,dimid2,dimid3/),varid)
       !
    ELSE IF(nbdim==3 .AND. typevar == 'double') THEN
       status = nf90_def_var(ncid,varname,nf90_double,     &
            (/dimid1,dimid2,dimid3/),varid)
       !
    ENDIF
    !
    status = nf90_enddef(ncid)
    IF(nbdim==4) status = nf90_put_var(ncid,varid,tabtemp4d)
    IF(nbdim==3) status = nf90_put_var(ncid,varid,tabtemp3d)
    !
    IF(ASSOCIATED( tabtemp3d ) ) DEALLOCATE( tabtemp3d )
    IF(ASSOCIATED( tabtemp4d ) ) DEALLOCATE( tabtemp4d )
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Write_Ncdf_var2d_Real_bis
  !
  !
  !
  !**************************************************************
  !   subroutine Write_Ncdf_var2d_real
  !**************************************************************
  SUBROUTINE Write_Ncdf_var2d_Real(varname,dimname,file,tabvar,typevar)
    !
	IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file,typevar
    CHARACTER(*), DIMENSION(2) :: dimname
    REAL*8, DIMENSION(:,:), POINTER :: tabvar
    !
    ! local variables
    !
    INTEGER :: dimid1,dimid2
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_WRITE,ncid)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_dimid(ncid,dimname(1), dimid1)
    status = nf90_inq_dimid(ncid,dimname(2), dimid2)
    status = nf90_inq_varid(ncid,varname,varid)
    status = nf90_redef(ncid)

    SELECT CASE(TRIM(typevar))
    CASE('double')
       status = nf90_def_var(ncid,varname,nf90_double,     &
            (/dimid1,dimid2/),varid)
    CASE('float')
       status = nf90_def_var(ncid,varname,nf90_float,     &
            (/dimid1,dimid2/),varid)
    END SELECT
    !
    status = nf90_enddef(ncid)
    status = nf90_put_var(ncid,varid,tabvar)
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Write_Ncdf_var2d_Real
  !
  !
  !
  !**************************************************************
  !   subroutine Write_Ncdf_var3d_real
  !**************************************************************
  SUBROUTINE Write_Ncdf_var3d_Real(varname,dimname,file,tabvar,typevar)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file,typevar
    CHARACTER(*),DIMENSION(3),INTENT(in) :: dimname
    REAL*8, DIMENSION(:,:,:), POINTER :: tabvar
    !
    ! local variables
    !
    INTEGER :: dimid1,dimid2,dimid3
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_WRITE,ncid)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_dimid(ncid,dimname(1), dimid1)
    status = nf90_inq_dimid(ncid,dimname(2), dimid2)
    status = nf90_inq_dimid(ncid,dimname(3), dimid3)
    status = nf90_inq_varid(ncid,varname,varid)
    status = nf90_redef(ncid)
    !
    SELECT CASE(TRIM(typevar))
    CASE('double')
       status = nf90_def_var(ncid,varname,nf90_double,     &
            (/dimid1,dimid2,dimid3/),varid)
    CASE('float')
       status = nf90_def_var(ncid,varname,nf90_float,     &
            (/dimid1,dimid2,dimid3/),varid)
    END SELECT
    !
    status = nf90_enddef(ncid)
    status = nf90_put_var(ncid,varid,tabvar)
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Write_Ncdf_var3d_Real
  !
  !
  !
  !**************************************************************
  !   subroutine Write_Ncdf_var4d_real
  !**************************************************************
  SUBROUTINE Write_Ncdf_var4d_Real(varname,dimname,file,tabvar,typevar)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file,typevar
    CHARACTER(*),DIMENSION(4),INTENT(in) :: dimname
    REAL*8, DIMENSION(:,:,:,:), POINTER :: tabvar
    !
    ! local variables
    !
    INTEGER :: dimid1,dimid2,dimid3,dimid4
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_WRITE,ncid)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_varid(ncid,varname,varid)
    !
    IF(status/=nf90_noerr) THEN
       !
       status = nf90_inq_dimid(ncid,dimname(1), dimid1)
       status = nf90_inq_dimid(ncid,dimname(2), dimid2)
       status = nf90_inq_dimid(ncid,dimname(3), dimid3)
       status = nf90_inq_dimid(ncid,dimname(4), dimid4)
       status = nf90_redef(ncid)
       !
       SELECT CASE(TRIM(typevar))
       CASE('double')
          status = nf90_def_var(ncid,varname,nf90_double,     &
               (/dimid1,dimid2,dimid3,dimid4/),varid)
       CASE('float')
          status = nf90_def_var(ncid,varname,nf90_float,     &
               (/dimid1,dimid2,dimid3,dimid4/),varid)
       END SELECT
       !
       status = nf90_enddef(ncid)
    ENDIF
    !
    status = nf90_put_var(ncid,varid,tabvar)
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Write_Ncdf_var4d_Real
  !
  !
  !
  !**************************************************************
  !   subroutine Write_Ncdf_var1d_Int
  !**************************************************************
  SUBROUTINE Write_Ncdf_var1d_Int(varname,dimname,file,tabvar)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file,dimname
    INTEGER, DIMENSION(:), POINTER :: tabvar
    !
    ! local variables
    !
    INTEGER :: dimid
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_WRITE,ncid)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_dimid(ncid,dimname, dimid)
    status = nf90_inq_varid(ncid,varname,varid)
    status = nf90_redef(ncid)
    status = nf90_def_var(ncid,varname,nf90_int,(/dimid/),varid)
    status = nf90_enddef(ncid)
    status = nf90_put_var(ncid,varid,tabvar)
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Write_Ncdf_var1d_Int
  !
  !
  !
  !**************************************************************
  !   subroutine Write_Ncdf_var2d_Int
  !**************************************************************
  SUBROUTINE Write_Ncdf_var2d_Int(varname,dimname,file,tabvar)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file
    CHARACTER(*), DIMENSION(2) :: dimname
    INTEGER, DIMENSION(:,:), POINTER :: tabvar
    !
    ! local variables
    !
    INTEGER :: dimid1,dimid2
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_WRITE,ncid)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_dimid(ncid,dimname(1), dimid1)
    status = nf90_inq_dimid(ncid,dimname(2), dimid2)
    status = nf90_inq_varid(ncid,varname,varid)
    status = nf90_redef(ncid)
    status = nf90_def_var(ncid,varname,nf90_int,     &
         (/dimid1,dimid2/),varid)
    status = nf90_enddef(ncid)
    status = nf90_put_var(ncid,varid,tabvar)
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Write_Ncdf_var2d_Int
  !
  !
  !
  !**************************************************************
  !   subroutine Write_Ncdf_var3d_Int
  !**************************************************************
  SUBROUTINE Write_Ncdf_var3d_Int(varname,dimname,file,tabvar)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file
    CHARACTER(*),DIMENSION(3),INTENT(in) :: dimname
    INTEGER, DIMENSION(:,:,:), POINTER :: tabvar
    !
    ! local variables
    !
    INTEGER :: dimid1,dimid2,dimid3
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_WRITE,ncid)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_dimid(ncid,dimname(1), dimid1)
    status = nf90_inq_dimid(ncid,dimname(2), dimid2)
    status = nf90_inq_dimid(ncid,dimname(3), dimid3)
    status = nf90_inq_varid(ncid,varname,varid)
    status = nf90_redef(ncid)
    status = nf90_def_var(ncid,varname,nf90_int,     &
         (/dimid1,dimid2,dimid3/),varid)
    status = nf90_enddef(ncid)
    status = nf90_put_var(ncid,varid,tabvar)
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Write_Ncdf_var3d_Int
  !
  !
  !
  !**************************************************************
  !   subroutine Write_Ncdf_var4d_Int
  !**************************************************************
  SUBROUTINE Write_Ncdf_var4d_Int(varname,dimname,file,tabvar)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file
    CHARACTER(*),DIMENSION(4),INTENT(in) :: dimname
    INTEGER, DIMENSION(:,:,:,:), POINTER :: tabvar
    !
    ! local variables
    !
    INTEGER :: dimid1,dimid2,dimid3,dimid4
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_WRITE,ncid)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_dimid(ncid,dimname(1), dimid1)
    status = nf90_inq_dimid(ncid,dimname(2), dimid2)
    status = nf90_inq_dimid(ncid,dimname(3), dimid3)
    status = nf90_inq_dimid(ncid,dimname(4), dimid4)
    status = nf90_inq_varid(ncid,varname,varid)
    status = nf90_redef(ncid)
    status = nf90_def_var(ncid,varname,nf90_int,     &
         (/dimid1,dimid2,dimid3,dimid4/),varid)
    status = nf90_enddef(ncid)
    status = nf90_put_var(ncid,varid,tabvar)
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Write_Ncdf_var4d_Int
  !
  !
  !
  !****************************************************************
  !               subroutine Read_Ncdf_var_t					*
  !						                                		*
  !  subroutine to read a variable in a given file for time t	*
  !						                                		*
  !     varname : name of variable to read		            	*
  !     file    : netcdf file name				                *
  !     tabvar  : values of the read variable	         		*
  !     time    : time corresponding to the values to read	    *
  !								                                *
  !****************************************************************
  !
  !
  !**************************************************************
  !          subroutine Read_Ncdf_var3d_real_t
  !**************************************************************
  SUBROUTINE Read_Ncdf_var3d_Real_t(varname,file,tabvar,time)
    !
    USE types
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file
    INTEGER,INTENT(in) :: time
    REAL*8, DIMENSION(:,:,:), POINTER :: tabvar
    !
    !local variables
    !
    INTEGER, DIMENSION(3) :: dimIDS
    INTEGER :: dim1,dim2
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_NOWRITE,ncid)
    !
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_varid(ncid,varname,varid)
    !
    status=nf90_inquire_variable(ncid,varid,dimids=dimIDS)
    status=nf90_inquire_dimension(ncid,dimIDS(1),len=dim1)
    status=nf90_inquire_dimension(ncid,dimIDS(2),len=dim2)
    !
    IF(.NOT. ASSOCIATED(tabvar)) THEN
       ALLOCATE(tabvar(dim1,dim2,1))
    ELSE
       IF( ANY(SHAPE(tabvar) /= (/dim1,dim2,1/)) ) THEN
          DEALLOCATE(tabvar)
          ALLOCATE(tabvar(dim1,dim2,1))
       ENDIF
    ENDIF

    status=nf90_get_var(ncid,varid,tabvar,start=(/1,1,time/))

    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to retrieve netcdf variable : ",TRIM(varname)
       STOP
    ENDIF
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Read_Ncdf_var3d_Real_t
  !
  !
  !
  !****************************************************************
  !              subroutine Write_Ncdf_var_t					*
  !		                                						*
  ! subroutine to write a variable in a given file for time t	*
  !		                                						*
  !     varname : name of variable to store		            	*
  !     dimname : name of dimensions of the given variable   	*
  !     file    : netcdf file name				                *
  !     tabvar  : values of the variable to write		    	*
  !     time    : time corresponding to the values to store	    *
  !								                                *
  !****************************************************************
  !
  !
  !**************************************************************
  !          subroutine Write_Ncdf_var3d_real_t
  !**************************************************************
  SUBROUTINE Write_Ncdf_var3d_Real_t(varname,dimname,file,tabvar,time,typevar)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file,typevar
    CHARACTER(*),DIMENSION(3),INTENT(in) :: dimname
    INTEGER :: time
    REAL*8, DIMENSION(:,:,:), POINTER :: tabvar
    !
    ! local variables
    !
    INTEGER :: dimid1,dimid2,dimid3
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_WRITE,ncid)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    IF(time==1) THEN
       !
       status = nf90_inq_dimid(ncid,dimname(1), dimid1)
       status = nf90_inq_dimid(ncid,dimname(2), dimid2)
       status = nf90_inq_dimid(ncid,dimname(3), dimid3)
       status = nf90_redef(ncid)

       !
       SELECT CASE(TRIM(typevar))
       CASE('double')
          status = nf90_def_var(ncid,varname,nf90_double,     &
               (/dimid1,dimid2,dimid3/),varid)
       CASE('float')
          status = nf90_def_var(ncid,varname,nf90_float,     &
               (/dimid1,dimid2,dimid3/),varid)
       END SELECT
       !
       status = nf90_enddef(ncid)

    ELSE
       status = nf90_inq_varid(ncid, varname, varid)
    ENDIF
    !
    status = nf90_put_var(ncid,varid,tabvar,start=(/1,1,time/))
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to store variable ",varname, &
	    " in file ",file
       STOP
    ENDIF
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Write_Ncdf_var3d_Real_t
  !
  !
  !
  !****************************************************************
  !               subroutine Read_Ncdf_var_t					*
  !			                                					*
  !  subroutine to read a variable in a given file for time t	*
  !					    	at level n	                        *
  !     varname : name of variable to read           			*
  !     file    : netcdf file name			                	*
  !     tabvar  : values of the read variable	        		*
  !     time    : time corresponding to the values to read   	*
  !     level   : level corresponding to the values to read	    *
  !								                                *
  !****************************************************************
  !
  !
  !**************************************************************
  !   subroutine Read_Ncdf_var4d_real_nt
  !**************************************************************
  SUBROUTINE Read_Ncdf_var4d_Real_nt(varname,file,tabvar,time,level)
    !
    USE types
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file
    INTEGER,INTENT(in) :: time,level
    REAL*8, DIMENSION(:,:,:,:), POINTER :: tabvar
    !
    !local variables
    !
    INTEGER, DIMENSION(4) :: dimIDS
    INTEGER :: dim1,dim2
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_NOWRITE,ncid)
    !
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_varid(ncid,varname,varid)
    !
    status=nf90_inquire_variable(ncid,varid,dimids=dimIDS)
    status=nf90_inquire_dimension(ncid,dimIDS(1),len=dim1)
    status=nf90_inquire_dimension(ncid,dimIDS(2),len=dim2)
    !
    IF(.NOT. ASSOCIATED(tabvar)) THEN
       ALLOCATE(tabvar(dim1,dim2,1,1))
    ELSE
       IF ((SIZE(tabvar,1) /= dim1) .OR. (SIZE(tabvar,2) /= dim2)) THEN
          DEALLOCATE(tabvar)
          ALLOCATE(tabvar(dim1,dim2,1,1))
       ENDIF
    ENDIF
    !
    status=nf90_get_var(ncid,varid,tabvar,start=(/1,1,level,time/),count=(/dim1,dim2,1,1/))
    !
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to retrieve netcdf variable : ",TRIM(varname)
       STOP
    ENDIF
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Read_Ncdf_var4d_Real_nt
  !
  !
  !
  !**************************************************************
  !   subroutine Read_Ncdf_var4d_real_t
  !**************************************************************
  SUBROUTINE Read_Ncdf_var4d_Real_t(varname,file,tabvar,time)
    !
    USE types
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file
    INTEGER,INTENT(in) :: time
    REAL*8, DIMENSION(:,:,:,:), POINTER :: tabvar
    !
    !local variables
    !
    INTEGER, DIMENSION(4) :: dimIDS
    INTEGER :: dim1,dim2,dim3
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_NOWRITE,ncid)
    !
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    status = nf90_inq_varid(ncid,varname,varid)
    !
    status=nf90_inquire_variable(ncid,varid,dimids=dimIDS)
    status=nf90_inquire_dimension(ncid,dimIDS(1),len=dim1)
    status=nf90_inquire_dimension(ncid,dimIDS(2),len=dim2)
    status=nf90_inquire_dimension(ncid,dimIDS(3),len=dim3)
    !
    IF(.NOT. ASSOCIATED(tabvar)) ALLOCATE(tabvar(dim1,dim2,dim3,1))
    status=nf90_get_var(ncid,varid,tabvar,start=(/1,1,1,time/))

    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to retrieve netcdf variable : ",TRIM(varname)
       STOP
    ENDIF
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Read_Ncdf_var4d_Real_t
  !
  !
  !
  !****************************************************************
  !               subroutine Write_Ncdf_var_t					*
  !				                                 				*
  !  subroutine to write a variable in a given file for time t	*
  !						    at level n	                        *
  !                                                             *
  !     varname : name of variable to store		            	*
  !     dimname : name of dimensions of the given variable   	*
  !     file    : netcdf file name				                *
  !     tabvar  : values of the variable to write			    *
  !     time    : time corresponding to the values to store	    *
  !     level   : level corresponding to the values to store	*
  !								                                *
  !****************************************************************
  !
  !
  !**************************************************************
  !   subroutine Write_Ncdf_var4d_real_t
  !**************************************************************
  SUBROUTINE Write_Ncdf_var4d_Real_t(varname,dimname,file,tabvar,time,typevar)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file,typevar
    CHARACTER(*),DIMENSION(4),INTENT(in) :: dimname
    INTEGER :: time,level
    REAL*8, DIMENSION(:,:,:,:), POINTER :: tabvar
    !
    ! local variables
    !
    INTEGER :: dimid1,dimid2,dimid3,dimid4
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_WRITE,ncid)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    IF(time==1) THEN
       !
       status = nf90_inq_dimid(ncid,dimname(1), dimid1)
       status = nf90_inq_dimid(ncid,dimname(2), dimid2)
       status = nf90_inq_dimid(ncid,dimname(3), dimid3)
       status = nf90_inq_dimid(ncid,dimname(4), dimid4)
       status = nf90_redef(ncid)
       !
       SELECT CASE(TRIM(typevar))
       CASE('double')
          status = nf90_def_var(ncid,TRIM(varname),nf90_double,     &
               (/dimid1,dimid2,dimid3,dimid4/),varid)
       CASE('float')
          status = nf90_def_var(ncid,TRIM(varname),nf90_float,     &
               (/dimid1,dimid2,dimid3,dimid4/),varid)
       END SELECT
       !
       status = nf90_enddef(ncid)

    ELSE
       status = nf90_inq_varid(ncid, varname, varid)
    ENDIF
    !
    status = nf90_put_var(ncid,varid,tabvar,start=(/1,1,1,time/))
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to store variable ",varname, &
	    " in file ",file
       STOP
    ENDIF
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Write_Ncdf_var4d_Real_t
  !
  !
  !
  !**************************************************************
  !   subroutine Write_Ncdf_var4d_real_nt
  !**************************************************************
  SUBROUTINE Write_Ncdf_var4d_Real_nt(varname,dimname,file,tabvar,time,level,typevar)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file,typevar
    CHARACTER(*),DIMENSION(4),INTENT(in) :: dimname
    INTEGER :: time,level
    REAL*8, DIMENSION(:,:,:,:), POINTER :: tabvar
    !
    ! local variables
    !
    INTEGER :: dimid1,dimid2,dimid3,dimid4
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_WRITE,ncid)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !
    IF(time==1.AND.level==1) THEN
       !
       status = nf90_inq_dimid(ncid,dimname(1), dimid1)
       status = nf90_inq_dimid(ncid,dimname(2), dimid2)
       status = nf90_inq_dimid(ncid,dimname(3), dimid3)
       status = nf90_inq_dimid(ncid,dimname(4), dimid4)
       status = nf90_redef(ncid)
       !
       SELECT CASE(TRIM(typevar))
       CASE('double')
          status = nf90_def_var(ncid,TRIM(varname),nf90_double,     &
               (/dimid1,dimid2,dimid3,dimid4/),varid)
       CASE('float')
          status = nf90_def_var(ncid,TRIM(varname),nf90_float,     &
               (/dimid1,dimid2,dimid3,dimid4/),varid)
       END SELECT
       !
       status = nf90_enddef(ncid)

    ELSE
       status = nf90_inq_varid(ncid, varname, varid)
    ENDIF
    !
    status = nf90_put_var(ncid,varid,tabvar,start=(/1,1,level,time/))
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to store variable ",varname, &
	    " in file ",file
       STOP
    ENDIF
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Write_Ncdf_var4d_Real_nt

  SUBROUTINE Write_Ncdf_var0d_Real(varname,file,tabvar,typevar)
    !
    IMPLICIT NONE
    !
    CHARACTER(*),INTENT(in) :: varname,file,typevar
    INTEGER :: time,level
    REAL*8 :: tabvar
    !
    ! local variables
    !
    INTEGER :: status,ncid
    INTEGER :: varid
    !
    status = nf90_open(file,NF90_WRITE,ncid)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",file
       STOP
    ENDIF
    !

    status = nf90_redef(ncid)
    !
    SELECT CASE(TRIM(typevar))
    CASE('double')
       status = nf90_def_var(ncid,TRIM(varname),nf90_double,     &
            varid=varid)
    CASE('float')
       status = nf90_def_var(ncid,TRIM(varname),nf90_float,     &
            varid=varid)
    END SELECT
    !
    status = nf90_enddef(ncid)

    !
    status = nf90_put_var(ncid,varid,tabvar)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to store variable ",varname, &
	    " in file ",file
       STOP
    ENDIF
    !
    status = nf90_close(ncid)
    !
  END SUBROUTINE Write_Ncdf_var0d_Real
  !
  !
  !****************************************************************
  !                 subroutine Read_Ncdf_VarName				*
  !			                                 					*
  !           subroutine to retrieve of all variables 			*
  !                  included in a given file					*
  !					                                			*
  !     filename    : netcdf file name	             			*
  !     tabvarname  : array containing various variables names	*
  !								                                *
  !****************************************************************
  !
  !**************************************************************
  !   subroutine Read_Ncdf_VarName
  !**************************************************************
  SUBROUTINE Read_Ncdf_VarName(filename,tabvarname)
    !
    CHARACTER(*),INTENT(in) :: filename
    CHARACTER*20,DIMENSION(:),POINTER :: tabvarname
    INTEGER :: nDimensions,nVariables
    INTEGER :: nAttributes,unlimitedDimId,i
    INTEGER :: ncid,status,dimid
    !
    status = nf90_open(filename,NF90_NOWRITE,ncid)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",filename
       STOP
    ENDIF
    !
    status = nf90_inquire(ncid,nDimensions,nVariables,nAttributes, &
         unlimitedDimId)
    !
    ALLOCATE(tabvarname(nVariables))
    !
    DO i=1,nVariables
       status = nf90_inquire_variable(ncid,i,tabvarname(i))
    END DO
    !
  END SUBROUTINE Read_Ncdf_Varname
  !
  !
  !
  !**************************************************************
  !   subroutine Copy_Ncdf_att
  !**************************************************************
  SUBROUTINE Copy_Ncdf_att_var(varname,filein,fileout)
    !
    CHARACTER(*),INTENT(in) :: filein,fileout
    CHARACTER(*),INTENT(in) :: varname
    INTEGER :: ncid_in,ncid_out,status,varid_in,varid_out
    !
    status = nf90_open(filein,NF90_NOWRITE,ncid_in)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open input netcdf file : ",filein
       STOP
    ENDIF
    !
    status = nf90_open(fileout,NF90_WRITE,ncid_out)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open output netcdf file : ",fileout
       STOP
    ENDIF
    !
    status = nf90_inq_varid(ncid_in,varname,varid_in)
    status = nf90_inq_varid(ncid_out,varname,varid_out)
    !
    status = nf90_redef(ncid_out)
    !
    status = nf90_copy_att(ncid_in,varid_in,'units',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'valid_min',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'valid_max',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'long_name',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'calendar',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'title',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'time_origin',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'positive',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'tstep_sec',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'nav_model',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'Minvalue=',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'Maxvalue=',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'short_name',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'online_operation',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'axis',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'interval_operation',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'interval_write',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'associate',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'actual_range',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'longitude',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'latitude',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'scale_factor',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'add_offset',ncid_out,varid_out)
    status = nf90_copy_att(ncid_in,varid_in,'missing_value',ncid_out,varid_out)
    !
    status = nf90_enddef(ncid_out)
    !
    status = nf90_close(ncid_in)
    status = nf90_close(ncid_out)
    !
  END SUBROUTINE Copy_Ncdf_att_var
  !
  !
  !
  !**************************************************************
  !       subroutine Copy_Ncdf_att
  !**************************************************************
  SUBROUTINE Copy_Ncdf_att_latlon(varname,filein,fileout,min,max)
    !
    CHARACTER(*),INTENT(in) :: filein,fileout
    CHARACTER(*),INTENT(in) :: varname
    REAL*8 :: min,max
    INTEGER :: ncid_in,ncid_out,status,varid_in,varid_out
    !
    status = nf90_open(filein,NF90_NOWRITE,ncid_in)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",filein
       STOP
    ENDIF
    !
    status = nf90_open(fileout,NF90_WRITE,ncid_out)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",fileout
       STOP
    ENDIF
    !
    status = nf90_inq_varid(ncid_in,varname,varid_in)
    status = nf90_inq_varid(ncid_out,varname,varid_out)
    !
    status = nf90_redef(ncid_out)
    !
    SELECT CASE (varname)
       !
    CASE('nav_lon')
       status = nf90_copy_att(ncid_in,varid_in,'units',ncid_out,varid_out)
       status = nf90_put_att(ncid_out,varid_out,'valid_min',REAL(min,4))
       status = nf90_put_att(ncid_out,varid_out,'valid_max',REAL(max,4))
       status = nf90_copy_att(ncid_in,varid_in,'long_name',ncid_out,varid_out)
       status = nf90_copy_att(ncid_in,varid_in,'nav_model',ncid_out,varid_out)
       status = nf90_copy_att(ncid_in,varid_in,'title',ncid_out,varid_out)
       !
    CASE('nav_lat')
       status = nf90_copy_att(ncid_in,varid_in,'units',ncid_out,varid_out)
       status = nf90_put_att(ncid_out,varid_out,'valid_min',REAL(min,4))
       status = nf90_put_att(ncid_out,varid_out,'valid_max',REAL(max,4))
       status = nf90_copy_att(ncid_in,varid_in,'long_name',ncid_out,varid_out)
       status = nf90_copy_att(ncid_in,varid_in,'nav_model',ncid_out,varid_out)
       status = nf90_copy_att(ncid_in,varid_in,'title',ncid_out,varid_out)
       !
    END SELECT
    !
    status = nf90_enddef(ncid_out)
    !
    status = nf90_close(ncid_in)
    status = nf90_close(ncid_out)
  END SUBROUTINE Copy_Ncdf_att_latlon
  !
  !
  !
  !**************************************************************
  !   function Get_NbDims
  !**************************************************************
  INTEGER FUNCTION Get_NbDims( varname , filename )
    !
    CHARACTER(*),INTENT(in) :: varname,filename
    INTEGER :: status,ncid,varid
    !
    status = nf90_open(TRIM(filename),NF90_NOWRITE,ncid)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",TRIM(filename)
       STOP
    ENDIF
    status = nf90_inq_varid(ncid,TRIM(varname),varid)
    status = nf90_inquire_variable(ncid, varid , ndims = Get_NbDims)
    !
    RETURN
    !
  END FUNCTION Get_NbDims
  !
  !
  !
  !**************************************************************
  !   function Get_NbDims
  !**************************************************************
  LOGICAL FUNCTION Dims_Existence( dimname , filename )
    !
    CHARACTER(*),INTENT(in) :: dimname,filename
    INTEGER :: status,ncid,dimid
    !
    status = nf90_open(TRIM(filename),NF90_NOWRITE,ncid)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",TRIM(filename)
       STOP
    ENDIF
    status = nf90_inq_dimid(ncid,dimname,dimid)
    !
    IF (status/=nf90_noerr) THEN
       Dims_Existence = .FALSE.
    ELSE
       Dims_Existence = .TRUE.
    ENDIF
    !
    RETURN
    !
  END FUNCTION Dims_Existence
  !
END MODULE io_netcdf


text:
      COMPLEX*16 FUNCTION javsubT0(Qs,Ts,mbt2,mtp2,mv2)
*------------------------------------------------------
      IMPLICIT NONE!
      REAL*8 Qs,Ts
      REAL*8 rmtp2,rmbt2,rmv2,zet2,eps,P,slamb2,x12,x22
      COMPLEX*16 mtp2,mbt2,mv2,C0,C02
      COMPLEX*16 YR,YP,YM,ieps,XSPENZ
      PARAMETER (zet2=1.6449340668482264364724151666460251D0)
      PARAMETER (eps=1D-30)
*
      rmtp2=DREAL(mtp2)
      rmbt2=DREAL(mbt2)
      rmv2 =DREAL(mv2)
      P=Qs+rmv2
      ieps=DCMPLX(0d0,eps)
      YR= (Qs+rmv2)/(Qs+rmtp2)
      YP= SQRT(ieps)
      YM=-SQRT(ieps)
*
      javsubT0=
     & -1d0/Ts*(
     & +DLOG(rmtp2*(-Qs)/(rmtp2+Qs)**2)*CDLOG((rmtp2-mv2)/(-Qs-mv2))
     &-CDLOG(rmtp2/(Ts-ieps)*(Ts/(rmtp2-rmv2)+1))*CDLOG(rmtp2/(Ts-ieps))
     & +CDLOG(rmtp2*(-P/rmtp2/(Ts-eps)-ieps))*CDLOG(rmtp2/(Ts-ieps))
     & -DLOG((rmtp2-rmv2)/rmtp2)*CDLOG(rmtp2/(Ts-ieps))
     & -1d0/2*CDLOG(-(Ts-ieps)/rmtp2)**2
     & +1d0/2*CDLOG(-(Ts-ieps)/rmtp2/(Ts/(rmtp2-rmv2)+1))**2
     & -(CDLOG(ieps/Ts-1d0-ieps/(rmtp2-rmv2))
     &  +CDLOG(ieps/Ts-rmv2/rmtp2+ieps/(rmtp2-rmv2))
     &  -CDLOG(P/rmtp2-ieps))*CDLOG((-Qs)/(Ts-ieps))
     & +CDLOG((P-ieps)/rmv2)*CDLOG((Qs-ieps)/(-Qs))
     & -1d0/2*CDLOG(rmv2/(Qs+rmv2-ieps))**2
     & -XSPENZ(rmv2/(Qs+rmv2-ieps))
     & +XSPENZ(1d0/(Ts-ieps)/(1/(rmtp2-rmv2)+1/Ts))
     & -XSPENZ(1d0/(rmv2/rmtp2-ieps/(rmtp2-rmv2)+ieps/Ts))
     & -XSPENZ(1d0/((rmtp2-rmv2)/rmtp2-ieps*(1/Ts-1/(rmtp2-rmv2))))
     &         )
*
     & +1d0/(rmtp2+Qs)*(
     &(-1d0/2*CDLOG(-1d0/((P-ieps)/(rmtp2-rmv2)))
     & -DLOG((-Qs)/(rmtp2+Qs)) )*CDLOG((-rmtp2+rmv2)/(P-ieps))
     & -CDLOG(-(1+Qs/rmtp2+ieps/(rmtp2-rmv2)))
     &               *DLOG(DABS((1+Qs/rmtp2)/(Qs/rmtp2)))
     & -CDLOG(-((rmv2+Qs)/rmtp2-ieps/(rmtp2-rmv2)))
     &               *DLOG(DABS((1+Qs/rmtp2)/(Qs/rmtp2)))
     & +CDLOG((P-ieps)/rmtp2*(1+Qs/rmtp2)-ieps)
     &               *DLOG(DABS((1+Qs/rmtp2)/(Qs/rmtp2)))
     & +XSPENZ((1+Qs/rmtp2)/(1+Qs/rmtp2+ieps/(rmtp2-rmv2)))
     & -XSPENZ((Qs/rmtp2)/(1+Qs/rmtp2+ieps/(rmtp2-rmv2)))
     & +XSPENZ((1+Qs/rmtp2)/((rmv2+Qs)/rmtp2-ieps/(rmtp2-rmv2)))
     & -XSPENZ((Qs/rmtp2)/((rmv2+Qs)/rmtp2-ieps/(rmtp2-rmv2)))
     & +XSPENZ(Qs/rmtp2/((1+Qs/rmtp2)-ieps))
     & -XSPENZ(DCMPLX(-1d0/((rmtp2+Qs)/(rmtp2-rmv2)-1),0d0))
     & +XSPENZ(1d0/((rmtp2-rmv2)/rmtp2+ieps/(rmtp2-rmv2)))
     &                 )
*
     & -(1d0/(rmtp2+Qs)+1d0/Ts)*(
     &+1d0/2*CDLOG((-rmtp2+mv2)/rmv2)
     &  *( CDLOG((-Qs-mv2)/rmv2)+2d0*CDLOG(-rmtp2/Qs*((Qs+mv2)/rmv2)**2)
     &    -CDLOG(-(-rmtp2+mv2)**2/rmv2**2))
     &-1d0/2*CDLOG((Qs+mv2)/rmv2)
     &     *(-CDLOG(-(Qs+mv2)*rmv2*rmv2/(Qs+rmv2)**3)
     &       +DLOG((+rmtp2-rmv2)/rmv2))
     &   -XSPENZ((Qs+mv2)/rmv2)
     &   +2d0*zet2
     &   +XSPENZ((-rmtp2+mv2)/rmv2)
     &   -XSPENZ((rmv2-rmtp2)/(Qs+rmtp2)/(YR-YP))
     &   -XSPENZ((rmv2-rmtp2)/(Qs+rmtp2)/(YR-YM))
     &                          )
*
      RETURN
      END

ext: f


text:
	SUBROUTINE DSGRP  ( igroup, iret )
C************************************************************************
C* DSGRP								*
C* 									*
C* This subroutine starts a new drawing element group.			*
C*									*
C* DSGRP  ( IGROUP, IRET )						*
C*									*
C* Input parameters:							*
C*	IGROUP		INTEGER		Group type			*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* S. Maxwell/GSC	 7/97						*
C************************************************************************
	INCLUDE 'FUNCCODE.PRM'
	INCLUDE 'ERROR.PRM'
C*
	INTEGER isend (3)
C------------------------------------------------------------------------
	isend (1) = 3
	isend (2) = CSGRP
	isend (3) = igroup
C
	CALL GPUT  ( isend, 3, iret )
	IF  ( iret .ne. NORMAL )  RETURN
C
	CALL GGET  ( iret , 1, ierr )
	IF  ( ierr .ne. NORMAL )  iret = ierr
C
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


program monte_carlo
      use types
      use mcrt_functions

      implicit none
      ! allocate variables and declare things
      real(dp), allocatable :: tau_list(:), total_scatter(:, :), mean_arr(:), angle_bins(:,:,:), position(:,:), ini_pos(:,:)
      real(dp) :: pos(3), nhat(3), old_pos(3)
      real(dp) :: r_max = 1.0, min_wt = 1e-3, albedo = 1.0_dp, p = 0.1
      real(dp) :: leng, r, tau, ran,  wt, exit_theta, exit_phi
      integer :: num_packets, len_tau, j, i, nbins = 200, u, phi_loc, theta_loc, k

      ! Parameters for the MC run
      num_packets = 1000000

      ! configure the tau that we want to see

      tau_list = [0.1, 1.0, 5.0, 10.0, 20.0, 50.0]
      len_tau = size(tau_list)

      ! allocating the memory for things
      allocate( total_scatter( len_tau, num_packets ) )
      allocate( mean_arr( len_tau ) )
      allocate( angle_bins( len_tau, nbins, nbins ) )
      allocate( position( num_packets, 3 ) )
      allocate( ini_pos( num_packets, 3 ) )

      ! begin program
      write(*,*) 'Program running!'

      do j = 1, size(tau_list)

            do i = 1, num_packets
                  wt = 1.0
                  ! r = 0.0_dp
                  tau = tau_list(j)

                  call blob_emit_bias(nhat)
                  call blob_initial_pos(pos, r_max) ! try to fix the separatrix issue and revert to old method?

                  leng = blob_initial_length(pos, nhat, r_max)
                  call mc_update(pos, nhat, leng)

                  ! force the particle to appear inside the sphere to trigger the MC simulation
                  ! this was done because of rounding errors
                  r = 0.0_dp
                  do while ((r < r_max) .and. (r >= 0.0_dp))
                        ! provided we are within the sphere
                        old_pos = pos
                        ! store the old position

                        ! generate the length
                        leng = mc_gen_L(tau, r_max)

                        ! create an nhat vector
                        if (r .ne. 0.0_dp) then ! convert into bool
                              call mc_emit(nhat)
                        end if

                        ! make it walk along the nhat vector a certain length leng
                        call mc_update(pos, nhat, leng)

                        ! if (i == 1) then
                        !       write(*,*) pos
                        ! end if

                        ! check if it will remain in the sphere
                        r = norm(pos)

                        if (r < r_max) then
                              ! calculate if we want to scatter further based on the albedo
                              ! implemented the russian roulette style giving packets an extra life before speaking to Kenny
                              ! he says it's not needed in this HW but I've included it because my simulations were run with it
                              if (ran .lt. albedo) then
                                    ! update the weight
                                    wt = albedo * wt
                                    ! let it scatter
                                    total_scatter(j,i) = total_scatter(j,i) + 1

                                    ! then do the killing of packets using the formula in the notes
                                    if (wt <= min_wt) then
                                    ! generate a random number
                                          call RANDOM_NUMBER(ran)
                                    ! either kill it or update the weight
                                          if (ran .lt. p) then
                                                wt = wt/p
                                          else
                                                continue
                                          end if
                                    end if
                              ! if we don't want to scatter then just let it continue without doing anything
                              else
                                    ! we might want to consider the fact that its then scattered straight through
                                    continue
                              end if

                        else if (r >= r_max) then
                              ! if it's outside the sphere then use the geometry to allow us to calculate the "true" length that it's travelled

                              leng = edge_length(old_pos, nhat, r_max)

                              ! now walk the particle to that location on the sphere, not to its final location
                              call mc_update(old_pos, nhat, leng)

                              ! calculate the exit angles in the far field based on the direction that the particle is going along
                              call exit_angles(nhat, exit_theta,  exit_phi)

                              ! generate a location in phi and theta as described in report
                              phi_loc = int( exit_phi/(pi)*nbins/2 + nbins/2 )
                              theta_loc = int( cos(exit_theta) * nbins/2 + nbins/2 +1)
                              ! bin the in the array
                              angle_bins(j, phi_loc, theta_loc) = angle_bins(j, phi_loc, theta_loc) + 1
                        end if
                  end do
            end do
            write(*,*) tau_list(j)
      end do

      open(newunit = u, file = "Q3/Data/angle_binning.txt", status = "replace")
      do k = 1, len_tau
              do j = 1, nbins
                  do i = 1, nbins
                        write(u,*) tau_list(k), ',', j,',',i , ',',angle_bins(k, j, i)
                  end do
            end do
      end do
end program


text:
C ***********************************************************
C
      DOUBLE PRECISION FUNCTION DONTR5(I,N,X,CONT,XSOL,HSOL)
C ----------------------------------------------------------
C     THIS FUNCTION CAN BE USED FOR CONINUOUS OUTPUT. IT PROVIDES AN
C     APPROXIMATION TO THE I-TH COMPONENT OF THE SOLUTION DERIVATIVE AT .
C     X. IT GIVES THE VALUE OF THE COLLOCATION POLYNOMIAL, DEFINED FOR
C     THE LAST SUCCESSFULLY COMPUTED STEP (BY DELAY5).
C ----------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER, PARAMETER :: DP=kind(1D0)
      REAL(kind=DP), dimension(1), intent(in) ::  CONT
C --- REQUIRED CONSTANTS
      COMMON /CONSTN/C1,C2,C1M1,C2M1,C1MC2

      N2=2*N
      N3=3*N

      S=(X-XSOL)/HSOL

	DONTR5=(CONT(I+N)+
     &   (S-C2M1)*(CONT(I+N2)+CONT(I+N3)*(S-C1M1))+
     &    S*(CONT(I+N2)+CONT(I+N3)*(2.*S-C1M1-C2M1)))/HSOL

      RETURN
      END
C
C     END OF FUNCTION DONTR5
C
C ***********************************************************


text:
!    ERGO: Efficient Robust Global Optimization
!    Copyright (C) 2007  Klaus Keller and Brian Tuttle
!
!    This program is free software; you can redistribute it and/or modify
!    it under the terms of the GNU General Public License as published by
!    the Free Software Foundation; either version 2 of the License, or
!    (at your option) any later version.
!
!    This program is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.
!
!    You should have received a copy of the GNU General Public License
!    along with this program; if not, write to the Free Software
!    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
!
!    Klaus Keller, kkeller@geosc.psu.edu
!    Brian Tuttle, btuttle@psu.edu
!    David McInerny, dmcinern@geosc.psu.edu
!    Alex Robinson, ajr225@psu.edu
!
!==========================================================================
!  This module provides an interface between the ERGO optimization package
!  and the numerical model that serves as an objective function.  Its
!  purpose is to centralize calls to the necessary model subroutines and
!  facilitate the incorporation of different models with the ERGO package.
!  The model itself may be defined here, or may be included as an established
!  module via a USE statement at the top with calls to its particular
!  subroutines in the appropriate places.
!==========================================================================
MODULE model

    USE global

    implicit none

    integer(i4b) :: nop
    real(DP), dimension(:), allocatable :: x_lower, x_upper
    real(DP) :: DE_pen_val

CONTAINS
!---------------------------------------------------------------------------
SUBROUTINE read_model_input(fid)
!  =========================================================================
! |  This subroutine is called from read_input() in module fileIO which     |
! |  opens the input file, reads the optimization section, and positions    |
! |  the pointer at the beginning of the model_specific section.  This      |
! |  routine continues reading the formatted data and assigns values to     |
! |  model variables.                                                       |
!  =========================================================================

    implicit none

! ID of a currently open file.
    integer(i4b), intent(IN) :: fid

    RETURN

END SUBROUTINE read_model_input
!---------------------------------------------------------------------------

!---------------------------------------------------------------------------
SUBROUTINE print_model_input()
!  ========================================================================
! |  This routine prints out the variables read from the input file and    |
! |  those computed in the model_init routine.  It is called from          |
! |  print_input() in ergo_fileIO.f90 and performs the model-specific      |
! |  part of that function.                                                |
!  ========================================================================

    implicit none

    RETURN

END SUBROUTINE print_model_input
!---------------------------------------------------------------------------

!---------------------------------------------------------------------------
SUBROUTINE model_write_arch(FID)
!  ========================================================================
! |  This routine writes the input variables to a file archive that can be |
! |  used as an input file to reproduce the calculation.  This is called   |
! |  from write_in_arch() in ergo_fileIO.f90.                              |
!  ========================================================================

    implicit none

! ID of a currently open file.
    integer(i4b), intent(IN) :: FID

    RETURN

END SUBROUTINE model_write_arch
!---------------------------------------------------------------------------

!---------------------------------------------------------------------------
SUBROUTINE model_size()
!  =========================================================================
! |  This subroutine conveys to the optimization package the sizes of the   |
! |  one-dimensional arrays for model parameters, penalty, and convergence. |
!  =========================================================================

    implicit none

! Default values for model sizes.
    n_param = 8
    n_param_pre = 1
    n_pen = 6
    n_pen_pre = 1
    n_conv = n_param
    n_conv_pre = n_param_pre

    RETURN

END SUBROUTINE model_size
!---------------------------------------------------------------------------

!---------------------------------------------------------------------------
SUBROUTINE model_init(eval_index, prereq_in)
!  ========================================================================
! |  This routine does any initialization that the model may require.  Any |
! |  prerequisite solution array may be included in the second argument.   |
! |  The output array is a list of parameter indices to be included by the |
! |  optimization program in evaluating the convergence of the solution.   |
! |  The default is to use all parameters for convergence evaluation.      |
!  ========================================================================

    implicit none

    integer(i4b), dimension(:), intent(OUT) :: eval_index
    real(DP), dimension(:), intent(IN), optional :: prereq_in
    integer :: i, astat

! Default value for pen_val in DE.
    DE_pen_val = 10000d0

! Initialize evaluation index with default values.
    eval_index = (/ (i, i=1, size(eval_index)) /)

    if (present(prereq_in)) then
        nop = n_param
        n_g = n_pen
    ! Allocate bounds.
        call model_allocate(nop)
    ! Initialize the model for the main optimization.

    else
        nop = n_param_pre
        n_g = n_pen_pre
    ! Allocate bounds.
        call model_allocate(nop)
    ! Initialize the model for the pre-calculation.

    end if

    x_lower = (/ 100d0,1000d0,1000d0,10d0,10d0,10d0,10d0,10d0 /)
    x_upper = (/10000d0,10000d0,10000d0,1000d0,1000d0,1000d0,1000d0,1000d0/)

    RETURN

END SUBROUTINE model_init
!---------------------------------------------------------------------------

!---------------------------------------------------------------------------
SUBROUTINE model_objf(x,f,g)
!  ========================================================================
! |  This is the model objective function.  The input is an array of       |
! |  decision parameters x.  The output is the single fitness value that   |
! |  rates the choice of parameters.  The constraint vector g provides a   |
! |  mechanism to penalize infeasible parameter choices.                   |
!  ========================================================================

    implicit none

    real(DP), dimension(nop), intent(IN) :: x
    real(DP), intent(OUT) :: f
    real(DP), dimension(n_g), intent(OUT) :: g
    integer :: nfe

! Initialize output.
    f = 0d0
    g = 0d0

! Default number of function evaluations is 1 for each call to model_objf.
    nfe = 1

! objective function
    f = x(1) + x(2) + x(3)

! constraints g<=0
    g(1) = -1.0d0 + 0.0025d0 * (x(4) + x(6))
    g(2) = -1.0d0 + 0.0025d0 * (x(5) + x(7) - x(4))
    g(3) = -1.0d0 + 0.01d0 * (x(8) - x(5))
    g(4) = -x(1) * x(6) + 833.33252d0 * x(4) + 100.0d0 * x(1) - 83333.333d0
    g(5) = -x(2) * x(7) + 1250.0d0 * x(5) + x(2) * x(4) - 1250.0d0 * x(4)
    g(6) = -x(3) * x(8) + 1250000.0d0 + x(3) * x(5) - 2500.0d0 * x(5)

    nfeval = nfeval + nfe

    RETURN

END SUBROUTINE model_objf
!---------------------------------------------------------------------------

!---------------------------------------------------------------------------
SUBROUTINE model_init_outfiles()
!  ========================================================================
! |  This routine provides an opportunity to initialize output files that  |
! |  will be updated by the model throughout the optimization process.     |
! |  File IDs for open files may be established here as well as header     |
! |  information for the output files.  It is called from init_outfiles()  |
! |  in ergo_fileIO.f90.                                                   |
!  ========================================================================

    implicit none

    RETURN

END SUBROUTINE model_init_outfiles
!---------------------------------------------------------------------------

!---------------------------------------------------------------------------
SUBROUTINE setup_test_objf(test_prereq_soln, test_solution)
!  ========================================================================
! |  Provide a known set of parameters to test the model objective         |
! |  function.                                                             |
!  ========================================================================

    implicit none

    real(DP), dimension(n_param_pre) :: test_prereq_soln
    real(DP), dimension(n_param) :: test_solution

    RETURN

END SUBROUTINE setup_test_objf
!---------------------------------------------------------------------------

!---------------------------------------------------------------------------
SUBROUTINE model_out(Solution, obf_out)
!  ========================================================================
! |  This output routine is called by the optimization program as it       |
! |  progresses toward convergence on a solution.  It allows the model to  |
! |  report the state of its various parameters given the latest best      |
! |  solution.  It is also called at the end of the optimization process.  |
!  ========================================================================

    implicit none

    real(DP), dimension(:), intent(IN) :: Solution
    real(DP), intent(OUT), OPTIONAL :: obf_out
    character(len=12) :: fname = "solution_out"
    real(DP), dimension(n_g) :: g
    real(DP) :: fout
    integer :: FID = 10
    integer :: istat, i

    call model_objf(Solution,fout,g)

    open(UNIT=FID, FILE=fname, STATUS='REPLACE', ACTION='WRITE', &
        ACCESS='SEQUENTIAL', IOSTAT=istat)

    write(FID,"(a,f11.4)") '# g10:  f(x) = ', fout
    write(FID,"(a)") '# i    x(i)'

    do i = 1,size(Solution)
        write(FID, "(i3,f11.4)") i, Solution(i)
    end do

    close(FID)

    if (present(obf_out)) obf_out = fout

    RETURN

END SUBROUTINE model_out
!---------------------------------------------------------------------------

!---------------------------------------------------------------------------
SUBROUTINE model_allocate(n_bounds)
!  ========================================================================
! |  Allocate upper and lower constraint arrays.                           |
!  ========================================================================

    implicit none

    integer, intent(IN) :: n_bounds
    integer :: astat

    astat = 0
    allocate(x_upper(n_bounds), STAT=astat)
    if (astat > 0) call error_msg(1,"x_upper")
    astat = 0
    allocate(x_lower(n_bounds), STAT=astat)
    if (astat > 0) call error_msg(1,"x_lower")

    RETURN

END SUBROUTINE model_allocate
!---------------------------------------------------------------------------

!---------------------------------------------------------------------------
SUBROUTINE model_finalize()
!  ========================================================================
! |  This routine cleans up memory allocated for model variables.          |
!  ========================================================================

    implicit none

    integer :: astat

! Deallocate upper and lower constraint arrays.
    astat = 0
    deallocate(x_upper, STAT=astat)
    if (astat > 0) call error_msg(2,"x_upper")
    astat = 0
    deallocate(x_lower, STAT=astat)
    if (astat > 0) call error_msg(2,"x_lower")

    RETURN

END SUBROUTINE model_finalize
!---------------------------------------------------------------------------
END MODULE model

max_stars_repo_name: fmyuan/clm-microbe

max_issues_repo_name: fmyuan/clm-microbe

max_forks_repo_name: email-clm/CLM-Microbe


text:
C
C	$Id: intzi.f,v 1.4 2008-07-27 00:17:22 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE INTZI (XX,YY,ZZ,LIN3,ITOP)
C
C FORCE STORAGE OF X, Y, AND Z INTO COMMON BLOCK
C
      COMMON /PWRZ2I/ X, Y, Z
      DATA IDUMX,IDUMY,IDUMZ /0, 0, 0/
      X = XX
      Y = YY
      Z = ZZ
      CALL INITZI (IDUMX,IDUMY,IDUMZ,LIN3,ITOP,1)
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
module erloff_call_stack_m
    use erloff_call_stack_entry_m, only: call_stack_entry_t
    use erloff_module_m, only: module_t
    use erloff_procedure_m, only: procedure_t
    use iso_varying_string, only: varying_string, operator(//)
    use strff, only: add_hanging_indentation, indent, join, NEWLINE

    implicit none
    private
    public :: call_stack_t

    type :: call_stack_t
        private
        type(call_stack_entry_t), allocatable :: entries(:)
    contains
        private
        procedure, public :: with_names_prepended
        procedure, public :: to_string
        procedure :: originated_from_module
        procedure :: originated_from_procedure
        generic, public :: operator(.originatedFrom.) => &
                originated_from_module, originated_from_procedure
        procedure :: includes_module
        procedure :: includes_procedure
        generic, public :: operator(.includes.) => &
                includes_module, includes_procedure
        procedure, public :: repr
    end type

    interface call_stack_t
        module procedure constructor
    end interface
contains
    pure function constructor(module_, procedure_) result(call_stack)
        type(module_t), intent(in) :: module_
        type(procedure_t), intent(in) :: procedure_
        type(call_stack_t) :: call_stack

        call_stack = internal_constructor([call_stack_entry_t(module_, procedure_)])
    end function

    pure function internal_constructor(entries) result(call_stack)
        type(call_stack_entry_t), intent(in) :: entries(:)
        type(call_stack_t) :: call_stack

        allocate(call_stack%entries, source = entries)
    end function

    pure function with_names_prepended(self, module_, procedure_) result(new_stack)
        class(call_stack_t), intent(in) :: self
        type(module_t), intent(in) :: module_
        type(procedure_t), intent(in) :: procedure_
        type(call_stack_t) :: new_stack

        new_stack = internal_constructor(&
                [call_stack_entry_t(module_, procedure_), self%entries])
    end function

    pure function to_string(self) result(string)
        class(call_stack_t), intent(in) :: self
        type(varying_string) :: string

        string = join(self%entries%to_string(), "->")
    end function

    pure function originated_from_module(self, module_)
        class(call_stack_t), intent(in) :: self
        type(module_t), intent(in) :: module_
        logical :: originated_from_module

        originated_from_module = self%entries(size(self%entries)).isFrom.module_
    end function

    pure function originated_from_procedure(self, procedure_)
        class(call_stack_t), intent(in) :: self
        type(procedure_t), intent(in) :: procedure_
        logical :: originated_from_procedure

        originated_from_procedure = self%entries(size(self%entries)).isFrom.procedure_
    end function

    pure function includes_module(self, module_)
        class(call_stack_t), intent(in) :: self
        type(module_t), intent(in) :: module_
        logical :: includes_module

        includes_module = any(self%entries.isFrom.module_)
    end function

    pure function includes_procedure(self, procedure_)
        class(call_stack_t), intent(in) :: self
        type(procedure_t), intent(in) :: procedure_
        logical :: includes_procedure

        includes_procedure = any(self%entries.isFrom.procedure_)
    end function

    pure function repr(self)
        class(call_stack_t), intent(in) :: self
        type(varying_string) :: repr

        repr = add_hanging_indentation( &
                "call_stack_t(" // NEWLINE &
                    // "entries = [" // NEWLINE &
                    // indent( &
                            join(self%entries%repr(), "," // NEWLINE), &
                            4) // NEWLINE // "]", &
                4) // NEWLINE // ")"
    end function
end module

max_stars_repo_name: everythingfunctional/erloff

max_issues_repo_name: everythingfunctional/erloff

max_forks_repo_name: everythingfunctional/erloff


text:
      subroutine doelem(elmcd)
      use maccbk
      use mackw
      use macttyp
      use macfile
      use macmisc
      implicit none
c
      integer*8 elmcd,ktcaloc
      character*(MAXSTR) token,wtoken*(*)
      integer slen,ival,ttype,slen2,ttype2,idx
      integer*4 hsrchz
      real*8 rval
      logical skipch
c
      call gettok(token,slen,ttype,rval,ival)
      if (.not. skipch(LCURL,token,slen,ttype,rval,ival) ) then
        if (ttype .eq. ttypID) then
c         call errmsg('doelem',
c    &              'missing '''//LCURL//''' is assumed',0,0)
        else
          call errmsg('doelem',
     &                ' syntax error:'//LCURL//' is expected.',0,16)
        endif
      endif
      goto 2000
C
c strat element definition
      entry doelm2(elmcd,wtoken,slen2,ttype2)
      slen=slen2
      ttype=ttype2
      token=wtoken(:slen)
 2000 if (skipch(COMMA,token,slen,ttype,rval,ival) ) go to 2000
c for debug
c      print *,'read element name ',ttype,ttypID,token(:slen)
c end debug
      if((token(:slen) .eq. RCURL) .or.
     &   (token(:slen) .eq. SEMIC))    then
         return
      else if (ttype .eq. ttypID) then
        idx =hsrchz(token(:slen))
        if (idtype(idx) .eq. icNULL) then
          idtype(idx)=int(elmcd)
          idval(idx)=ktcaloc(kytbl(kwMAX,elmcd)+1)
c          if(elmcd .eq. icMARK)then
c            write(*,*)'doelem ',idx,idval(idx)
c          endif
          call rdkwdl(idx)
c          if(elmcd .eq. icMARK)then
c            write(*,*)'doelem-read ',rlist(idval(idx)+1)
c          endif

c for debug
c      print *,'doelem:',idtype(idx),elmcd,pname(kytbl(0,elmcd))
c    &         ,kytbl(kwMAX,elmcd)
c end debug
        else if(idtype(idx) .eq. elmcd) then
c         call errmsg('doelm',
c    &              'element: '//token(:slen)//' is redefined',0,0)
          call rdkwdl(idx)
        else if(idtype(idx) .lt. icMXEL) then
          call errmsg('doelem',
     &      'syntax error: '//token(:slen)//
     $         ' is already used as '//
     &      pname(kytbl(0,idtype(idx)))//'.'
     &        ,0,16)
        else if (idtype(idx) .gt. icLINE) then
          call errmsg('doelem',
     &        token(:slen)//'is used as a element name.'
     &        ,0,0)
          call errmsg('doelem',
     &        'syntax error:you can not use reserved words as a name.'
     &        ,0,16)
        else
          call errmsg('doelem',
     &      'syntax error: '//token(:slen)//' is already used as'//
     &      ' anohter element name.'
     &        ,0,16)
        endif
      else
        call errmsg('doelm',
     &       'syntax error:doelm{<line name>=(line definition)}',0,16)
      endif
c
      call gettok(token,slen,ttype,rval,ival)
      go to 2000
c
      end

max_stars_repo_name: noboruatkek/SAD

max_issues_repo_name: noboruatkek/SAD

max_forks_repo_name: noboruatkek/SAD


! { dg-do run }
!
! Check the fix for PR67779, in which array sections passed in the
! recursive calls to 'quicksort' had an incorrect offset.
!
! Contributed by Arjen Markus  <arjen.markus895@gmail.com>
!
! NOTE: This is the version of the testcase in comment #16 (from Thomas Koenig)
!
module myclass_def
    implicit none

    type, abstract :: myclass
    contains
        procedure(assign_object), deferred        :: copy
        procedure(one_lower_than_two), deferred   :: lower
        procedure(print_object), deferred         :: print
        procedure, nopass                         :: quicksort  ! without nopass, it does not work
    end type myclass

    abstract interface
        subroutine assign_object( left, right )
            import                        :: myclass
            class(myclass), intent(inout) :: left
            class(myclass), intent(in)    :: right
        end subroutine assign_object
    end interface

    abstract interface
        logical function one_lower_than_two( op1, op2 )
            import                     :: myclass
            class(myclass), intent(in) :: op1, op2
        end function one_lower_than_two
    end interface

    abstract interface
        subroutine print_object( obj )
            import                     :: myclass
            class(myclass), intent(in) :: obj
        end subroutine print_object
    end interface

    !
    ! Type containing a real
    !

    type, extends(myclass) :: mysortable
        integer :: value
    contains
        procedure :: copy    => copy_sortable
        procedure :: lower   => lower_sortable
        procedure :: print   => print_sortable
    end type mysortable

contains
!
! Generic part
!
recursive subroutine quicksort( array )
    class(myclass), dimension(:) :: array

    class(myclass), allocatable :: v, tmp
    integer                     :: i, j

    integer :: k

    i = 1
    j = size(array)

    allocate( v,   source = array(1) )
    allocate( tmp, source = array(1) )

    call v%copy( array((j+i)/2) ) ! Use the middle element

    do
        do while ( array(i)%lower(v) )
            i = i + 1
        enddo
        do while ( v%lower(array(j)) )
            j = j - 1
        enddo

        if ( i <= j ) then
            call tmp%copy( array(i) )
            call array(i)%copy( array(j) )
            call array(j)%copy( tmp )
            i        = i + 1
            j        = j - 1
        endif

        if ( i > j ) then
            exit
        endif
    enddo

    if ( 1 < j ) then
        call quicksort( array(1:j) ) ! Problem here
    endif

    if ( i < size(array) ) then
        call quicksort( array(i:) )  ! ....and here
    endif
end subroutine quicksort

!
! Specific part
!
subroutine copy_sortable( left, right )
    class(mysortable), intent(inout) :: left
    class(myclass), intent(in)       :: right

    select type (right)
        type is (mysortable)
            select type (left)
                type is (mysortable)
                    left = right
            end select
    end select
end subroutine copy_sortable

logical function lower_sortable( op1, op2 )
    class(mysortable), intent(in) :: op1
    class(myclass),    intent(in) :: op2

    select type (op2)
        type is (mysortable)
            lower_sortable = op1%value < op2%value
    end select
end function lower_sortable

subroutine print_sortable( obj )
    class(mysortable), intent(in) :: obj

    write(*,'(G0," ")', advance="no") obj%value
end subroutine print_sortable

end module myclass_def

! test program
program test_quicksort
    use myclass_def

    implicit none

    type(mysortable), dimension(20) :: array
    real, dimension(20) :: values

    call random_number(values)

    array%value = int (1000000 * values)

! It would be pretty perverse if this failed!
    if (check (array)) STOP 1

    call quicksort( array )

! Check the the array is correctly ordered
    if (.not.check (array)) STOP 2
contains
     logical function check (arg)
         type(mysortable), dimension(:) :: arg
         integer                        :: s
         s = size (arg, 1)
         check = all (arg(2 : s)%value .ge. arg(1 : s - 1)%value)
     end function check
end program test_quicksort


text:
!#########################################################
!#########################################################
!#########################################################
!#########################################################
subroutine rho_ana(x,d,dx,ncell)
  use amr_parameters
  use hydro_parameters
  use poisson_parameters
  implicit none
  integer ::ncell                         ! Number of cells
  real(dp)::dx                            ! Cell size
  real(dp),dimension(1:nvector)       ::d ! Density
  real(dp),dimension(1:nvector,1:ndim)::x ! Cell center position.
  !================================================================
  ! This routine generates analytical Poisson source term.
  ! Positions are in user units:
  ! x(i,1:3) are in [0,boxlen]**ndim.
  ! d(i) is the density field in user units.
  !================================================================
  integer::i
  real(dp)::emass,xmass,ymass,zmass,rr,rx,ry,rz,rcut
  real(dp)::v200,c,spin,fgas,dnfw,dmax,fcut

  emass=2.*boxlen*0.5d0**nlevelmax
  xmass=boxlen/2.0
  ymass=boxlen/2.0
  zmass=boxlen/2.0
  rcut =boxlen/2.0*0.75

  v200=gravity_params(1)
  c   =gravity_params(2)
  spin=gravity_params(3)
  fgas=gravity_params(4)
  dmax=1.0/emass/(1.0+emass)**2

  do i=1,ncell
     rx=x(i,1)-xmass
     ry=x(i,2)-ymass
     rz=x(i,3)-zmass
     rr=sqrt(rx**2+ry**2+rz**2)
     fcut=exp(-(rr/rcut)**50)
     dnfw=1d0/rr/(1d0+rr)**2
     dnfw=MIN(dmax,dnfw)*fcut
     dnfw=MAX(dnfw,1d-15)
     d(i)=(1d0-fgas)*dnfw
  end do

end subroutine rho_ana

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rknop/amuse