Varfit Subroutine

SUBROUTINE Varfit &
!
(varModel, modelselected)

IMPLICIT NONE

!arguments with intent (in):
INTEGER (KIND = short), INTENT(IN) :: varModel

!arguments with intent(out):
INTEGER (KIND = short), INTENT(OUT) :: modelselected !!if automodel is passed as varModel, the optimal model is returned

!local declarations:
REAL (KIND = float) :: ami1, ama1, ami2, ama2, ami3, ama3 
REAL (KIND = float) :: ami4, ama4, ami5, ama5
REAL (KIND = float) :: gmax, gmin, dmin, dmax
REAL (KIND = float) :: fac, fape
INTEGER :: i,j, k
INTEGER (KIND = short) :: nmodels, model
REAL (KIND = float) :: sillmodel (3), rangemodel (3), nuggetmodel (3), of (3), iterations (3)
REAL (KIND = float) :: minimumOF !!minimum value of objective function
!-----------------------------------end of declarations-----------------------!

!!find minimum and maximum values
dmin = MINVAL ( dist, MASK = dist > 0.0 )
dmax = MAXVAL ( dist)
gmin = MINVAL ( empVariogram, MASK = empVariogram > 0.0   )
gmax = MAXVAL ( empVariogram  )

ipara = 2

IF ( ndir > 1 ) ipara = 4
IF ( nvar == 1 ) ipara = ipara - 1
IF ( npep == 1 ) ipara = ipara + 1

!set number of sevivariogram models to fit
IF (varmodel == AUTOFIT) THEN
    nmodels = 2
ELSE
    nmodels = 1
END IF

DO k = 1, nmodels
    
    !reset variables
    ieva = 0 
    nugget = 0.0
    range = 0.0
    partialSill = 0.0
    anisotropyAngle = 0.0
    anisotropyRatio = 1.0
    var = expVar
    
    IF (nmodels == 1) THEN
        model = varmodel
    ELSE
        model = k
    END IF
    
        !     ALCANCE MAYOR
              AMI1=0.1*DMIN
              AMA1=1.5*DMAX
        !     ANGULO DE ANISOTROPIA
              AMI2=0.0
              AMA2=PI
        !     RATIO DE ANISOTROPIA
              AMI3=1.0
              AMA3=20.0
        !     MESETA
              AMI4=0.1*VAR
              AMA4=2.0*VAR
        !     EFECTO DE PEPITA
              AMI5=0.0
              AMA5=VAR
              IF (model == 4) THEN
                FAC=2.0*GMAX/DMAX
                FAPE=2.0*GMIN
              END IF
        !        WRITE (6,*)'PARAMETERS TO ESTIMATE: ',IPARA
                parRangeMin(1)=AMI1
                parRangeMax(1)=AMA1
                IF (IPARA.EQ.2) THEN
                  IF (NPEP.EQ.1) THEN
                    parRangeMin(2)=AMI5
                    parRangeMax(2)=AMA5
                  ELSE
                    parRangeMin(2)=AMI4
                    parRangeMax(2)=AMA4
                    IF (varModel.EQ.4) THEN
                      parRangeMin(2)=0.00001
                      parRangeMax(2)=FAC
                    END IF
                  END IF
                ELSE IF (IPARA.EQ.3) THEN
                  IF (NPEP.EQ.1) THEN
                    parRangeMin(2)=AMI4
                    parRangeMax(2)=AMA4
                    parRangeMin(3)=AMI5
                    parRangeMax(3)=AMA5
                    IF (model.EQ.4) THEN
                      parRangeMin(2)=0.00001
                      parRangeMax(2)=FAC
                      parRangeMin(3)=0.0
                      parRangeMax(3)=FAPE
                    END IF
                  ELSE
                    parRangeMin(2)=AMI2
                    parRangeMax(2)=AMA2
                    parRangeMin(3)=AMI3
                    parRangeMax(3)=AMA3
                  END IF
                ELSE IF (IPARA.EQ.4) THEN
                  parRangeMin(2)=AMI2
                  parRangeMax(2)=AMA2
                  parRangeMin(3)=AMI3
                  parRangeMax(3)=AMA3
                  IF (NPEP.EQ.1) THEN
                    parRangeMin(4)=AMI5
                    parRangeMax(4)=AMA5
                  ELSE
                    parRangeMin(4)=AMI4
                    parRangeMax(4)=AMA4
                    IF (model.EQ.4) THEN
                      parRangeMin(4)=0.00001
                      parRangeMax(4)=FAC
                    END IF
                  END IF
                ELSE
                  parRangeMin(2)=AMI2
                  parRangeMax(2)=AMA2
                  parRangeMin(3)=AMI3
                  parRangeMax(3)=AMA3
                  parRangeMin(4)=AMI4
                  parRangeMax(4)=AMA4
                  parRangeMin(5)=AMI5
                  parRangeMax(5)=AMA5
                  IF (model.EQ.4) THEN
                    parRangeMin(4)=0.00001
                    parRangeMax(4)=FAC
                    parRangeMin(5)=0.0
                    parRangeMax(5)=FAPE
                  END IF
                END IF

        !!     VARIANCE = SILL + NUGGET

              IF (model.EQ.4) THEN
                 parRangeMin(1)=0.00001
                 parRangeMax(1)=1.9
              END IF   

              CALL Simplex (model)
              
              sillmodel (k) = partialSill
              rangemodel (k) = range
              nuggetmodel (k) = nugget
              of (k) = hfina
              iterations (k) = ieva
END DO



IF (nmodels > 1) THEN !set the best model
    minimumOF = 10E6
     DO i = 1, nmodels
          IF ( of (i) < minimumOF ) THEN
              minimumOF = of (i)
              modelselected = i
              partialSill = sillmodel (i)
              range = rangemodel (i)
              nugget = nuggetmodel (i)
              ieva = iterations (i)
              
          END IF
     END DO
ELSE
    modelselected = varModel
END IF

        
IF (nmodels > 1) THEN
    hfina = minimumOF
END IF

RETURN
END SUBROUTINE Varfit