Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.

Lahey Fortran TO Intel Fortran

Javier_C_3
Beginner
1,149 Views

Hello,

I'm finishing my PhD in Geotechincal Engineering, in the Polytechnic University of Catalonia.

I've implemented Fortran Code within Plaxis environtment (User Soil Defined Model). To create the necessaries dynamic libraries .dll, I've been using Intel Fortran Compiler (32/64) and LF95.

Intel Fortran (32/64) compiles the principal part of the code, and Lahey compiles an additional part of the code that let me use Plaxis Interface introducing input parameters. 

I've using next make orders (which I've adapted from the original ones) for compiling the code and create the .dll:

IFORT32 (ifcomp.bat)
echo !DEC$ ATTRIBUTES DLLExport, StdCall, reference :: User_Mod
ifort @comp_rel  UsrMod.f90 UsrAddDF.f90 usrlib.f90 okmsgivf.f90 /dll /exe:example_ivf

IFORT64 (ifcomp64.bat)
REM compile 64-bit DLL: make sure you're in a 64-bit development environment
copy impivf impexp
ifort /static  UsrMod.f90 usrlib.f90 okmsgivf.f90 /dll /exe:example_ivf64
 

LF95 (lfcomp.bat)
echo       DLL_Export User_Mod
lf95 -c usrlib.f90
lib /out:lfusrlib.lib usrlib.ob
lf95 -c -ml winapi okmsglah.f90
lf95 -win -dll  UsrMod.f90 Usr_Add.f90 okmsglah.obj -lib lfusrlib.lib -ml bd -out:example_lah.dll


I've been able to compile the code of the constitutive code that I've developed in PLAXIS, introduce the model parameters, calculate and analyze the results. But now, LF95 trial license has expired (>30 days), and Lahey doesn't offer me any Student Version at a reasonable price.

 

 

That's why I'd like to ask you how to change LF95 make orders to Intel Fortran. I've tried all it was in my hand, but I'm not an expert developer (can manage well numerical models, but get lost with compiler orders -except basic ones-).


Thank you very much.
Best regards,

 

0 Kudos
11 Replies
mecej4
Honored Contributor III
1,149 Views

The requirements for the Fortran DLL are on page 174 of the Plaxis documentation http://www.plaxis.nl/files/files/2D-3-Material-Models.pdf .

You need to add the DEC$ directive given on the same page to your Fortran source code, and use the compiler command listed on the same page.

 Mentioning the Lahey compiler in this context is an unnecessary distraction.

0 Kudos
Javier_C_3
Beginner
1,149 Views

Thak you for your answer.

The structure of the compiler orders that I've written in the previous message is the same as the one you could find in Plaxis USRMOD example in their web (they use both compilers):"example_ivf.dll" and "example_ivf64.dll" are created with Intel Compiler and "example_lah.dll" with Lahey Compiler. "example_ivf.dll" and "example_ivf64.dll" are correct.

I've tried many times to change the directives of lfcomp.bat (the one that creates "example_lah.dll") to Intel directives, but I failed (Plaxis didn't work). That's why I asked.

echo       DLL_Export User_Mod        ==>  echo !DEC$ ATTRIBUTES DLLExport, StdCall, reference :: User_Mod
lf95 -c usrlib.f90                                ==>   ¿?¿? (userlib.f90 have no DEC of DLL_Export directive)
lib /out:lfusrlib.lib usrlib.obj                 ==> ¿?¿?  
lf95 -c -ml winapi okmsglah.f90           ==> ¿?¿?  (there is an equivalent code with DEC directives, 
okmsgivf.f90)
lf95 -win -dll  UsrMod.f90 Usr_Add.f90 okmsglah.obj -lib lfusrlib.lib -ml bd -out:example_lah.dll   ==> If I write any of the following orders, it doesn'st work:

ifort /winapp UsrMod.f90  dfusrlib.lib /dll   (not work)

ifort /winapp UsrMod.f90  usrlib.lib /dll   (not work)

ifort /winapp UsrMod.f90  usrlib.90 /dll   (not work)

ifort /winapp UsrMod.f90 UsrAddDF.f90 usrlib.lib /dll   (not work)

ifort /winapp UsrMod.f90 UsrAddDF.f90 okmsgivf.f90 usrlib.lib /dll   (not work)

ifort /winapp UsrMod.f90 UsrAddDF.f90 okmsgivf.f90 usrlib.f90 /dll   (not work)

ifort /winapp UsrMod.f90 UsrAddDF.f90 usrlib.f90 /dll /exe:example_lah  (not work)
 
ifort @comp_rel UsrMod.f90 UsrAddDF.f90 usrlib.f90 /dll /exe:example_lah  (not work)

 

Comp_rel include: /real_size:64
/fpconstant
/optimize:4
/static
/assume:noaccuracy_sensitive
/assume:byterecl
/assume:buffered_io
/notraceback
/debug:none
/check:none
/heap-arrays:1
/fpp
/recursive

ifort @comp_rel UsrMod.f90 UTsrAddDF.f90 okmsgivf.f90 usrlib.lib /dll /exe:example_lah  (not work)

Etc...

 

Thank You

0 Kudos
mecej4
Honored Contributor III
1,149 Views

I don't understand what the problem is.

Apparently, you have Intel Fortran and you are able to produce with it DLLs that Plaxis can use correctly? And, you used to have Lahey F95 and you were able to use it (until the trial license expired) to produce DLLs which also worked correctly with Plaxis? If so, what is holding you back from creating IFort-compatible versions of the sources that you used with LF95?

What do you mean by "not work" and "it doesn't work"? Failure to compile, failure to link, rejection of DLL by Plaxis? In any case, you should show us the details of the attempted commands and the error messages that you received.

Note that you have not shown any source code. Plaxis is a big commercial package and I am not going to be able to do anything that requires having Plaxis or Delphi. If you post the files UsrMod.f90, UTsrAddDF.f90, okmsgivf.f90 and usrlib.f90, we can help with building an IFort DLL with them.

 

 

0 Kudos
Javier_C_3
Beginner
1,149 Views

The problem is that I don't know how to "translate" the make file (and some lines in the code) from LF95 to Intel.

I'm sure that it's possible to work only with Intel, but the files that Plaxis provides as an example of User Defined Soil Model (http://kb.plaxis.nl/models/plaxis-user-defined-soil-models) uses both compilers, and, as I'm not a professional developer, I've followed this structure.

Example Plaxis file uses:

- INTEL to create two .dll:  example_ivf.dll (compiled with IA-32) and example_ivf64.dll (compiled with IA-64)
- LF95: to create one .dll: example_lah.dll

I changed ".for" => ".f90", and adapted makefiles. With that, I implemented a very simple model, and after finishing compilation, I put this 3 files in Plaxis folder and all worked well.

In order to use LF95 or Intel Compiler, it's necessary to change one line of the code of  USRMOD.f90 (see below).

After that, I started the implementation of an elastoplastic model which already exist in Plaxis (hardening soil). I'm copying below the code of this one.

When I finish its validation (making comparations with Plaxis Hardening Soil Model), I will implement my model code (which tries to improve the elastic part of that one).

I don't have any problem to copy here the code, but it has more than 15.000 lines, so I'm going to copy the main files and some subroutines (if you need more information, please, tell me).

Thank you very much.


The code I'm attaching works perfectly once compiled with INTEL and LF95 (example_ivf.dll with IA-32, example_ivf64.dll with IA-64, and example_lah.dll with LF95).

It uses the libraries provided by Plaxis, and Plaxis structure for subrotuines. As well, it has some errors from a physical point of view which I've already correct but not implemented yet (in any case, they're not important for our propose here).

I compiled this code with next make files (which I executed from the respectives compiler consoles):

IFCOMP.BAT: 
echo !DEC$ ATTRIBUTES DLLExport, StdCall, reference :: User_Mod
ifort @comp_rel  UsrMod.f90 UsrAddDF.f90 usrlib.f90 okmsgivf.f90 /dll /exe:example_ivf

COMP_REL:
/real_size:64
/fpconstant
/optimize:4
/static
/assume:noaccuracy_sensitive
/assume:byterecl
/assume:buffered_io
/notraceback
/debug:none
/check:none
/heap-arrays:1
/fpp
/recursive

 

IFCOMP64.BAT:
REM compile 64-bit DLL: make sure you're in a 64-bit development environment
copy impivf impexp
ifort /static  UsrMod.f90 usrlib.f90 okmsgivf.f90 /dll /exe:example_ivf64

 

IMPIVF:
: User_Mod

LFCOMP.BAT:
echo       DLL_Export User_Mod
lf95 -c usrlib.f90
lib /out:lfusrlib.lib usrlib.obj
lf95 -c -ml winapi okmsglah.f90
lf95 -win -dll  UsrMod.f90 Usr_Add.f90 okmsglah.obj -lib lfusrlib.lib -ml bd -out:example_lah.dll

CODE:

USRMOD.F90 (TO USE WITH INTEL OR LF95)
_____________________________________________________________________________________________________

      Subroutine User_Mod ( IDTask, iMod, IsUndr, &
                           iStep, iTer, iEl, Int, &
                           X, Y, Z, &
                           Time0, dTime, &
                           Props, Sig0, Swp0, StVar0, &
                           dEps, D, BulkW, &
                           Sig, Swp, StVar, ipl, &
                           nStat, NonSym, iStrsDep, iTimeDep,iTang, &
                           iPrjDir, iPrjLen, iAbort )

!
      Implicit Double Precision (A-H, O-Z)
!
      Dimension Props(*), Sig0(*), StVar0(*), dEps(*), D(6,6), &
                Sig(*), StVar(*), iPrjDir(*)
      Character*255 PrjDir, Dbg_Name
      Data iounit / 0 /
      Save iounit
      
      ! La siguiente linea debe utilizarse para la compilación del DLL 
      !NEXT TWO LINES MUST BE ACTIVATE OR DESACTIVATE DEPENDING ON THE COMPILER USED:
      !DEC$ ATTRIBUTES DLLExport, StdCall, Reference :: User_Mod  ! activate with Intel IA-32/64
      !DLL_Export User_Mod                                        !activate with LF95
 
      ! Posibilidad de debugging 
      If (iounit.Eq.0) Then
        PrjDir=' '
        Do i=1,iPrjLen
          PrjDir(i:i) = Char(iPrjDir(i))
        End Do

        Dbg_Name=PrjDir(:iPrjLen)//'udsmex'
        nErr=0
    1   Continue
          Open(Unit= 1, File= Dbg_Name, iostat=ios)
          If (ios.Eq.0) Close(Unit=1, Status='delete', iostat=ios)

          If (ios.Ne.0) Then
            !
            ! in case of error try ...udsmex1 or udsmex2 or ..
            !
            nErr=nErr+1
            dbg_name=PrjDir(:iPrjLen)//'udsmex'//char(48+nErr)

            If (nErr.Lt.10) Goto 1
          End If

        Open(Unit= 1, File=Dbg_Name, blocksize=4096)

        Write(1,*)'File 1 opened'
        iounit = 1
!        Call WriVec(-1,'Props',Props,14)
      End If
!      Call WriIvl(-1, 'iounit',iounit )
!      Call WriIvl(-1, 'IDTask',IDTask )
      Select Case (iMod)
        Case (1)   
          Call MyModel1( IDTask, iMod, IsUndr, iStep, iTer, iEl, Int, &
                        X, Y, Z, Time0, dTime, &
                        Props, Sig0, Swp0, StVar0, &
                        dEps, D, BulkW, Sig, Swp, StVar, ipl, &
                        nStat, NonSym, iStrsDep, iTimeDep, iTang, &
                        iAbort )
!        Case (2)
!          Call MyModel3( IDTask, ....
        Case Default
          Write(1,*) 'invalid model number in UsrMod', iMod
          Write(1,*) 'IDTask: ',IDTask
          Stop 'invalid model number in UsrMod'
          iAbort=1
          Return
      End Select ! iMod
      If (IDTask .EQ. 5 .AND. iel+int .EQ. 2) Then
        Write(1,*)'nStat   : ',nStat
        Write(1,*)'NonSym  : ',NonSym
        Write(1,*)'StrsDep : ',iStrsDep
        Write(1,*)'TimeDep : ',iTimeDep
        Write(1,*)'Tangent : ',iTang
      End If
      If (IDTask .EQ. 3 .AND. iel+int .EQ. 2) Then
        Write(1,*)'IDTask: ', IDTask,'iStep: ', iStep, 'iTer :', iTer
        Call Flush(1)
      End If
!      Call WriIvl(-1, 'IDTask end',IDTask )
      Return
      End ! User_Mod

      include 'Matriz6.f90'      
      include 'Matriz3.f90'        
      include 'mymodel1.f90'
      include 'HS_plast.f90'
      include 'MC_plast.f90'      
      include 'TRAC_plast.f90'
      include 'CAP_plast.f90'      
      include 'HS_CAP_plast.f90'
      include 'MC_CAP_plast.f90'
      include 'HS_TRAC_plast.f90'
      include 'MC_TRAC_plast.f90'      
 _____________________________________________________________________________________________________ 

 

USR_ADD.F90  (TO USE WITH LF95)
_____________________________________________________________________________________________________

! Subroutines in this file:

! Exported:
!  Subroutine GetModelCount( nMod )
!  Subroutine GetModelName ( iMod , ModelName )
!  Subroutine GetParamCount( iMod , nParam )
!  Subroutine GetParamName ( iMod , iParam, ParamName )
!  Subroutine GetParamUnit ( iMod , iParam, Units )
!  Subroutine GetStateVarCount( iMod , nVar )
!  Subroutine GetStateVarName ( iMod , iVar, Name )
!  Subroutine GetStateVarUnit ( iMod , iVar, Unit )
!
! Local:
!  Subroutine GetParamAndUnit( iMod , iParam, ParamName, Units )
!  Subroutine GetStateVarNameAndUnit( iMod , iVar, Name, Unit )
!
      Subroutine GetModelCount(nMod)
      !
      ! Return the maximum model number (iMod) in this DLL
      !
      Integer (Kind=4) nMod
      Dll_Export GetModelCount

      nMod = 1 ! Maximum model number (iMod) in current DLL

      Return
      End ! GetModelCount

      Subroutine GetModelName( iMod , ModelName )
      !
      ! Return the name of the different models
      !
      Integer  iMod
      Character (Len= * ) ModelName
      Character (Len=255) tName
      DLL_Export GetModelName

      Select Case (iMod)
        Case (1)
          tName = 'HS'
        Case Default
          tName = 'not in DLL'
      End Select
      LT = Len_Trim(tName)
      ModelName= Char(lt) // tName(1:Lt)
      Return
      End ! GetModelName

      Subroutine GetParamCount( iMod , nParam )
      !
      ! Return the number of parameters of the different models
      !
      DLL_Export GetParamCount
      Select Case (iMod)
        Case ( 1 )
          nParam = 14
        Case Default
          nParam = 0
      End Select
      Return
      End ! GetParamCount

      Subroutine GetParamName( iMod , iParam, ParamName )
      !
      ! Return the parameters name of the different models
      !
      Character (Len=255) ParamName, Units
      DLL_Export GetParamName
      Call GetParamAndUnit(iMod,iParam,ParamName,Units)
      Return
      End

      Subroutine GetParamUnit( iMod , iParam, Units )
      !
      ! Return the units of the different parameters of the different models
      !
      Character (Len=255) ParamName, Units
      DLL_Export GetParamUnit
      Call GetParamAndUnit(iMod,iParam,ParamName,Units)
      Return
      End

      Subroutine GetParamAndUnit( iMod , iParam, ParamName, Units )
      !
      ! Return the parameters name and units of the different models
      !
      ! Units: use F for force unit
      !            L for length unit
      !            T for time unit
      !
      Character (Len=255) ParamName, Units, tName
      Select Case (iMod)
        Case (1)
          ! ModName = 'EL'
          Select Case (iParam)
            Case (1)
              ParamName = 'G'          ; Units     = 'F/L^2#'
            Case (2)
              ParamName = 'nu'         ; Units     = '-'
            Case (3)
              ParamName = 'E50_ref'    ; Units     = 'F/L^2#'
            Case (4)
              ParamName = 'Eoed_ref'   ; Units     = 'F/L^2#'
            Case (5)
              ParamName = 'Eur_ref'    ; Units     = 'F/L^2#'
            Case (6)
              ParamName = 'm'          ; Units     = '-'
            Case (7)
              ParamName = 'c'          ; Units     = 'F/L^2#'
            Case (8)
              ParamName = 'phi'        ; Units     = 'deg'
            Case (9)
              ParamName = 'psi'        ; Units     = 'deg'
            Case (10)
              ParamName = 'nu_ur'      ; Units     = '-'
            Case (11)
              ParamName = 'p_ref'      ; Units     = 'F/L^2#'            
            Case (12)
              ParamName = 'K0_nc'      ; Units     = '-'
            Case (13)
              ParamName = 'R_f'        ; Units     = '-'              
            Case (14)
              ParamName = 'sitma_tracc'; Units     = 'F/L^2#'              
            Case Default
              ParamName = '???'     ; Units     = '???'
          End Select
        Case Default
          ! model not in DLL
          ParamName = ' N/A '        ; Units     = ' N/A '
      End Select
      tName    = ParamName
      LT       = Len_Trim(tName)
      ParamName= Char(lt) // tName(1:Lt)

      tName = Units
      LT    = Len_Trim(tName)
      Units = Char(lt) // tName(1:Lt)

      Return
      End ! GetParamAndUnit


      Subroutine GetStateVarCount( iMod , nVar )
      !
      ! Return the number of state variables of the different models
      !

      Dll_export GetStateVarCount

      	Select Case (iMod)
        Case (1)
          nVar = 2
        Case Default
          nVar = 0
      End Select

      Return
      End

      Subroutine GetStateVarName( iMod , iVar, Name )
      !
      ! Return the name of the different state variables
      ! of the different models
      !

      Character (Len=255) Name, Unit

      Dll_export GetStateVarName

        Call GetStateVarNameAndUnit( iMod , iVar, Name, Unit )

      End

      Subroutine GetStateVarUnit( iMod , iVar, Unit )
      !
      ! Return the units of the different state variables of the different models
      !
      Character (Len=255) Name, Unit

      Dll_export GetStateVarUnit

        Call GetStateVarNameAndUnit( iMod , iVar, Name, Unit )

      End

      Subroutine GetStateVarNameAndUnit( iMod , iVar, Name, Unit )
      !
      ! Return the name and unit of the different state variables of the different models
      !
      Character (Len=255) Name, Unit, tName
      	Select Case (iMod)
        Case (1)
          Select Case (iVar)
            Case (1)
              Name = 'gamma_p'         ; Unit = '-'
            Case (2)
              Name = 'p_p'             ; Unit = 'F/L^2#'             
            Case Default
              Name='N/A'              ; Unit = '?'
          End Select
        Case Default
          Name='N/A'                  ; Unit = '?'
      End Select
      tName    = Name
      LT       = Len_Trim(tName)
      Name     = Char(lt) // tName(1:Lt)

      tName    = Unit
      LT       = Len_Trim(tName)
      Unit     = Char(lt) // tName(1:Lt)
      Return
      End
_____________________________________________________________________________________________________
 
 

USRADDDF.F90  (TO USE WITH INTEL)
_____________________________________________________________________________________________________

! Subroutines in this file:
! Exported:
!  Subroutine GetModelCount( nMod )
!  Subroutine GetModelName ( iMod , ModelName )
!  Subroutine GetParamCount( iMod , nParam )
!  Subroutine GetParamName ( iMod , iParam, ParamName )
!  Subroutine GetParamUnit ( iMod , iParam, Units )
!  Subroutine GetStateVarCount( iMod , nVar )
!  Subroutine GetStateVarName ( iMod , iVar, Name )
!  Subroutine GetStateVarUnit ( iMod , iVar, Unit )
!
! Local:
!  Subroutine GetParamAndUnit( iMod , iParam, ParamName, Units )
!  Subroutine GetStateVarNameAndUnit( iMod , iVar, Name, Unit )


      Subroutine GetModelCount(nMod)
      !
      ! Return the maximum model number (iMod) in this DLL
      !
      Integer (Kind=4) nMod
!DEC$ ATTRIBUTES DLLExport, StdCall, reference ::  GetModelCount

      nMod = 1 ! Maximum model number (iMod) in current DLL

      Return
      End ! GetModelCount

      Subroutine GetModelName( iMod , ModelName )
      !
      ! Return the name of the different models
      !
      Integer  iMod
      Character (Len= * ) ModelName
      Character (Len=255) tName
!DEC$ ATTRIBUTES DLLExport, StdCall, reference ::  GetModelName

      Select Case (iMod)
        Case (1)
          tName = 'HS'
        Case Default
          tName = 'not in DLL'
      End Select
      LT = Len_Trim(tName)
      ModelName= Char(lt) // tName(1:Lt)
      Return
      End ! GetModelName

      Subroutine GetParamCount( iMod , nParam )
      !
      ! Return the number of parameters of the different models
      !
!DEC$ ATTRIBUTES DLLExport, StdCall, reference ::  GetParamCount
      Select Case (iMod)
        Case ( 1 )
          nParam = 14
        Case Default
          nParam = 0
      End Select
      Return
      End ! GetParamCount

      Subroutine GetParamName( iMod , iParam, ParamName )
      !
      ! Return the parameters name of the different models
      !
      Character (Len=255) ParamName, Units
!DEC$ ATTRIBUTES DLLExport, StdCall, reference ::  GetParamName
      Call GetParamAndUnit(iMod,iParam,ParamName,Units)
      Return
      End

      Subroutine GetParamUnit( iMod , iParam, Units )
      !
      ! Return the units of the different parameters of the different models
      !
      Character (Len=255) ParamName, Units
!DEC$ ATTRIBUTES DLLExport, StdCall, reference ::  GetParamUnit
      Call GetParamAndUnit(iMod,iParam,ParamName,Units)
      Return
      End

      Subroutine GetParamAndUnit( iMod , iParam, ParamName, Units )
      !
      ! Return the parameters name and units of the different models
      !
      ! Units: use F for force unit
      !            L for length unit
      !            T for time unit
      !
      Character (Len=255) ParamName, Units, tName
      Select Case (iMod)
        Case (1)
          ! ModName = 'HS'
          Select Case (iParam)
            Case (1)
              ParamName = 'G'          ; Units     = 'F/L^2#'
            Case (2)
              ParamName = 'nu'         ; Units     = '-'
            Case (3)
              ParamName = 'E50_ref'    ; Units     = 'F/L^2#'
            Case (4)
              ParamName = 'Eoed_ref'   ; Units     = 'F/L^2#'
            Case (5)
              ParamName = 'Eur_ref'    ; Units     = 'F/L^2#'
            Case (6)
              ParamName = 'm'          ; Units     = '-'
            Case (7)
              ParamName = 'c'          ; Units     = 'F/L^2#'
            Case (8)
              ParamName = 'phi'        ; Units     = 'deg'
            Case (9)
              ParamName = 'psi'        ; Units     = 'deg'
            Case (10)
              ParamName = 'nu_ur'      ; Units     = '-'
            Case (11)
              ParamName = 'p_ref'      ; Units     = 'F/L^2#'            
            Case (12)
              ParamName = 'K0_nc'      ; Units     = '-'
            Case (13)
              ParamName = 'R_f'        ; Units     = '-'     
            Case (14)
              ParamName = 'sitma_tracc'; Units     = 'F/L^2#'                 
            Case Default
              ParamName = '???'   ; Units     = '???'
          End Select
        Case Default
          ! model not in DLL
          ParamName = ' N/A '     ; Units     = ' N/A '
      End Select
      tName    = ParamName
      LT       = Len_Trim(tName)
      ParamName= Char(lt) // tName(1:Lt)

      tName = Units
      LT    = Len_Trim(tName)
      Units = Char(lt) // tName(1:Lt)

      Return
      End ! GetParamAndUnit

      Subroutine GetStateVarCount( iMod , nVar )
      !
      ! Return the number of state variables of the different models
      !

!DEC$ ATTRIBUTES dllexport, stdcall, reference :: GetStateVarCount

      	Select Case (iMod)
        Case (1)
          nVar = 2
        Case Default
          nVar = 0
      End Select

      Return
      End

      Subroutine GetStateVarName( iMod , iVar, Name )
      !
      ! Return the name of the different state variables of the different models
      !
      Character (Len=255) Name, Unit

!DEC$ ATTRIBUTES dllexport, stdcall, reference :: GetStateVarName

        Call GetStateVarNameAndUnit( iMod , iVar, Name, Unit )

      End

      Subroutine GetStateVarUnit( iMod , iVar, Unit )
      !
      ! Return the units of the different state variables of the different models
      !
      Character (Len=255) Name, Unit

!DEC$ ATTRIBUTES dllexport, stdcall, reference :: GetStateVarUnit

        Call GetStateVarNameAndUnit( iMod , iVar, Name, Unit )
      End

      Subroutine GetStateVarNameAndUnit( iMod , iVar, Name, Unit )
      !
      ! Return the name and unit of the different state variables of the different models
      !
      Character (Len=255) Name, Unit, tName
      	Select Case (iMod)
        Case (1)
          Select Case (iVar)
            Case (1)
              Name = 'gamma_p'         ; Unit = '-'
            Case (2)
              Name = 'p_p'             ; Unit = 'F/L^2#'              
            Case Default
              Name='N/A'              ; Unit = '?'
          End Select
        Case Default
          Name='N/A'                  ; Unit = '?'
      End Select
      tName    = Name
      LT       = Len_Trim(tName)
      Name     = Char(lt) // tName(1:Lt)

      tName    = Unit
      LT       = Len_Trim(tName)
      Unit     = Char(lt) // tName(1:Lt)
      Return
      End
_____________________________________________________________________________________________________

 

OKMSGIVF.F90  (TO USE WITH INTEL)
_____________________________________________________________________________________________________
      Subroutine OK_MessageBox(t)
      use dfwin
      use dfwinty
      integer (kind=Int_Ptr_Kind()) hWnd

      character*(*) t

      character(len=256) :: msg, title

!  Display a messagebox with an OK button
!  Note that all strings must be null terminated for C's sake

      msg    = Trim(t)  // char(0)
      title = 'UDSM_Example' // char(0)

      hWnd = 0

      iret = MessageBox( hWnd, &
                        msg, &
                        title, &
                        MB_OK )

      Return
      End
_____________________________________________________________________________________________________
 
 
OKMSGLAH.F90  (TO USE WITH LF95)
_____________________________________________________________________________________________________
      Subroutine OK_MessageBox(t)
      dll_import MessageBoxA

      integer :: MessageBoxA
      Integer is3DF

!  Some (not all) Message box flags. The C header file winuser.h
!  contains a complete list of flags. Winuser.h is not provided
!  in the LF95 distribution.

      integer, parameter :: MB_OK                   = 0

      integer            :: iret
      character(len=*)   :: t
      character(len=256) :: msg, title

!  Display a messagebox with an OK button
!  Note that all strings must be null terminated for C's sake

      msg  = Trim(t)  // char(0)
      title = 'UDSM_example' // char(0)

      iret = MessageBoxA(val(0),&
                        val(pointer(msg)), &
                        val(pointer(title)), &
                        val(MB_OK))

      Return
      End
_____________________________________________________________________________________________________
 
 
USRLIB.F90  (SUBROUTINES; LIBRARY)
_____________________________________________________________________________________________________
Subroutine MZEROR(R,K)
!C
!C***********************************************************************
!C
!C     Function: To make a real array R with dimension K to zero
!C
!C***********************************************************************
!C
      Implicit Double Precision (A-H,O-Z)
      Dimension R(*)

      Do J=1,K
        R(J) = 0.0D0
      End Do

      Return
      End


      Subroutine MZEROI(I,K)
!C
!C***********************************************************************
!C
!C     Function: To make an integre array I with Dimension K to zero
!C
!C***********************************************************************
!C
      Dimension I(*)

      Do J=1,K
        I(J)=0
      End Do

      Return
      End

      Subroutine SETRVAL(R,K,V)
!C
!C***********************************************************************
!C
!C     Function: To fill a real array R with Dimension K with value V
!C
!C***********************************************************************
!C
      Implicit Double Precision (A-H,O-Z)
      Dimension R(*)

      Do J=1,K
        R(J)=V
      End Do

      Return
      End

      Subroutine SETIVAL(I,K,IV)
!C
!C***********************************************************************
!C
!C     Function: To fill an integer array I with Dimension K with value IV
!C
!C***********************************************************************
!C
      Implicit Double Precision (A-H,O-Z)
      Dimension I(*)

      Do J=1,K
        I(J)=IV
      End Do

      Return
      End

      Subroutine COPYIVEC(I1,I2,K)
!C
!C***********************************************************************
!C
!C     Function: To copy an integer array I1 with Dimension K to I2
!C
!C***********************************************************************
!C
      Implicit Double Precision (A-H,O-Z)
      Dimension I1(*),I2(*)

      Do  J=1,K
        I2(J)=I1(J)
      End Do

      Return
      End

      Subroutine COPYRVEC(R1,R2,K)
!C
!C***********************************************************************
!C
!C     Function: To copy a Double array R1 with Dimension K to R2
!C
!C***********************************************************************
!C
      Implicit Double Precision (A-H,O-Z)
      Dimension R1(*),R2(*)

      Do J=1,K
        R2(J)=R1(J)
      End Do

      Return
      End


      Logical Function IS0ARR(A,N)
!C
!C***********************************************************************
!C    Function :  To check whether a real array contains only zero values.
!C                When an array contains only zero's is might not need to be
!C                written to the XXX file.
!C                exit Function when first non-zero value occured or when
!C                all elements are checked and are zero.
!C
!C    Input:  A : array to be checked
!C            N : number of elements in array that should be checked
!C
!!C    Output : .TRUE.  when all elements are 0
!C             .FALSE. when at least one element is not zero
!C
!C    Called by :  Subroutine TOBXX
!C
!C***********************************************************************
!C
      Implicit Double Precision (A-H,O-Z)
      Dimension A(*)
      Is0Arr=.False.
      Do I=1,N
        If ( A(I) .Ne. 0 ) Return
      End Do
      Is0Arr=.True.
      Return
      End

      Logical Function IS0IARR(IARR,N)
!C
!C***********************************************************************
!C    Function :  To check whether a integer array contains only zero values.
!C                Similar to IS0ARR
!C
!C    Input:  IARR : array to be checked
!C            N    : number of elements in array that should be checked
!C
!C    Output : .TRUE.  when all elements are 0
!C             .FALSE. when at least one element is not zero
!C
!C    Called by :  Subroutine TOBXX
!C
!C***********************************************************************
!C
      Implicit Double Precision (A-H,O-Z)
      Dimension IARR(*)

      Is0IArr=.False.
      Do I=1,N
        If ( IARR(I) .Ne. 0 ) Return
      End Do
      Is0IArr=.True.
      Return
      End
!C***********************************************************************
      Subroutine MulVec(V,F,K)
!C***********************************************************************
!C
!C     Function: To multiply a real vector V with dimension K by F
!!C
!C***********************************************************************
!C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION V(*)

      Do J=1,K
        V(J)=F*V(J)
      End Do

      Return
      End     ! Subroutine Mulvec
!C***********************************************************************
      Subroutine MatVec(xMat,IM,Vec,N,VecR)
!!C***********************************************************************
!C
!C     Calculate VecR = xMat*Vec
!!C
!C I   xMat  : (Square) Matrix (IM,*)
!C I   Vec   : Vector
!C I   N     : Number of rows/colums
!C O   VecR  : Resulting vector
!C
!C***********************************************************************
      Implicit Double Precision (A-H,O-Z)
      Dimension xMat(IM,*),Vec(*),VecR(*)
!C***********************************************************************
      Do I=1,N
        X=0
        Do J=1,N
          X=X+xMat(I,J)*Vec(J)
        End Do
        VecR(I)=X
      End Do
      Return
      End    ! Subroutine MatVec

!C***********************************************************************
      Subroutine AddVec(Vec1,Vec2,R1,R2,N,VecR)
!C***********************************************************************
!C
!C     Calculate VecR() = R1*Vec1()+R2*Vec2()
!C
!C I   Vec1,
!C I   Vec2  : Vectors
!C I   R1,R2 : Multipliers
!C I   N     : Number of rows
!C O   VecR  : Resulting vector
!C
!C***********************************************************************
      Implicit Double Precision (A-H,O-Z)
      Dimension Vec1(*),Vec2(*),VecR(*)
!C***********************************************************************
      Do I=1,N
        X=R1*Vec1(I)+R2*Vec2(I)
        VecR(I)=X
      End Do
      Return
      End    ! Subroutine AddVec
!C
!C***********************************************************************
      Subroutine DInProd(A,B,N,RES)
!C***********************************************************************
!C
!C     Returns the Inproduct of two vectors
!C
!C I   A,B  : Two vectors
!C I   N    : Used length of vectors
!C***********************************************************************
      Implicit Double Precision (A-H,O-Z)
      Dimension A(*),B(*)
!C***********************************************************************

      X = 0
      Do I=1,N
        X = X + A(I)*B(I)
      End Do
      RES = X
      Return
      End     ! Function DInProd
!C
!C***********************************************************************
      Subroutine MatMat(xMat1,Id1,xMat2,Id2,nR1,nC2,nC1,xMatR,IdR)
!C***********************************************************************
!!C
!C     Calculate xMatR = xMat1*xMat2
!C
!C I   xMat1 : Matrix (Id1,*)
!C I   xMat2 : Matrix (Id2,*)
!C I   nR1   : Number of rows in resulting matrix    (= No rows in xMat1)
!C I   nC2   : Number of columns in resulting matrix (= No cols in xMat2)
!C I   nC1   : Number of columns in matrix xMat1
!C             = Number  rows    in matrix xMat2
!C O   xMatR : Resulting matrix (IdR,*)
!C
!C***********************************************************************
      Implicit Double Precision (A-H,O-Z)
      Dimension xMat1(Id1,*),xMat2(Id2,*),xMatR(IdR,*)
!C**********************************************************************

      Do I=1,nR1
        Do J=1,nC2
          X=0
          Do K=1,nC1
            X=X+xMat1(I,K)*xMat2(K,J)
          End Do
          xMatR(I,J)=X
        End Do
      End Do

      Return
      End     ! Subroutine MatMat

!C***********************************************************************
      Subroutine MatMatSq(n, xMat1, xMat2, xMatR)
!C***********************************************************************
!C
!C     Calculate xMatR = xMat1*xMat2 for square matrices, size n
!C
!C I   n     : Dimension of matrices
!C I   xMat1 : Matrix (n,*)
!C I   xMat2 : Matrix (n,*)
!C O   xMatR : Resulting matrix (n,*)
!C
!C***********************************************************************
      Implicit Double Precision (A-H,O-Z)
      Dimension xMat1(n,*),xMat2(n,*),xMatR(n,*)
!C**********************************************************************

      Do I=1,n
        Do J=1,n
          X=0
          Do K=1,n
            X=X+xMat1(I,K)*xMat2(K,J)
          End Do
          xMatR(I,J)=X
        End Do
      End Do

      Return
      End     ! Subroutine MatMatSq

!C***********************************************************************
      Subroutine WriVal ( io, C , V )
!C***********************************************************************
!C
!!C Write (Double) value to file unit io (when io>0)
!C
!C***********************************************************************
!C
      Implicit Double Precision (A-H,O-Z)
      Character C*(*)

      If (io.Le.0) Return

      Write(io,*) C,V
    1 Format( A,3x, 1x,1p,e12.5)
      Return
      End
!C***********************************************************************
      Subroutine WriIVl ( io, C , I )
!C***********************************************************************
!C
!C Write (integer) value to file unit io (when io>0)
!C
!C***********************************************************************
      Implicit Double Precision (A-H,O-Z)
      Character C*(*)

      If (io.Le.0) Return

      Write(io,*) C,I
    1 Format( A,3x, 1x,I6)
      Return
      End
!C***********************************************************************
      Subroutine WriIVc ( io, C , iV , n )
!C***********************************************************************
!C
!C Write (integer) vector to file unit io (when io>0)
!C
!C***********************************************************************
      Character C*(*)
      Dimension iV(*)

      If (io.Le.0) Return

      Write(io,*) C
      Write(io,1) (iv(i),i=1,n)
    1 Format( ( 2(3x,5i4) ) )
      Return
      End
!C***********************************************************************
      Subroutine WriVec ( io, C , V , n )
!C***********************************************************************
!C
!C Write (Double) vector to file unit io (when io>0)
!C 6 values per line
!C***********************************************************************
      Implicit Double Precision (A-H,O-Z)
      Character C*(*)
      Dimension V(*)

      If (io.Le.0) Return

      If (Len_Trim(C).Le.6) Then
        Write(io,2) C,( V(i),i=1,n)
      Else
        Write(io,*) C
        Write(io,1) ( V(i),i=1,n)
      End If
    1 Format( ( 2(1x, 3(1x,1p,e10.3) ) ) )
    2 Format( A, ( T7, 2(1x, 3(1x,1p,e10.3) ) ) )
      Return
      End
!C***********************************************************************
      Subroutine WriVec5( io, C , V , n )
!C***********************************************************************
!C
!C Write (Double) vector to file unit io (when io>0)
!C 5 values per line
!C***********************************************************************
      Implicit Double Precision (A-H,O-Z)
      Character C*(*)
      Dimension V(*)

      If (io.Le.0) Return

      Write(io,*) C
      Write(io,1) ( V(i),i=1,n)
    1 Format( 5(1x,1p,e12.5) )
      Return
      End
!C***********************************************************************
      Subroutine WriMat ( io, C , V , nd, nr, nc )
!C***********************************************************************
!C
!C Write (Double) matrix to file unit io (when io>0)
!C 6 values per line
!C***********************************************************************
      Implicit Double Precision (A-H,O-Z)
      Character C*(*)
      Dimension V(nd,*)

      If (io.Le.0) Return

      Write(io,*) C
      Do j=1,nr
        Write(io,1) j,( V(j,i),i=1,nc)
      End Do
    1 Format(i4, (  T7,2(1x, 3(1x,1p,e10.3) ) ) )
      Return
      End
!C***********************************************************************

      Subroutine MatInvPiv(Aorig,B,N)
      Implicit Double Precision (A-H,O-Z)
      Dimension Aorig(n,*), B(n,*),A(:,:)
      Allocatable :: A
      Allocate ( A(n,n) ) ! No error checking !!
      Call CopyRVec(AOrig, A, n*n )
      Call MZeroR(B,n*n)
      Do i=1,n
        B(i,i) = 1d0
      End Do
      Do I=1,n
        T=A(I,I)
        iPiv=i
        Do j=i+1,n
          If ( Abs(A(j,i)) .Gt. Abs(A(iPiv,i))  ) iPiv=j
        End Do
        If (iPiv.Ne.i) Then
          Do j=1,n
            x         = A( i  ,j)
            A( i  ,j) = A(iPiv,j)
            A(iPiv,j) = x
            x         = B( i  ,j)
            B( i  ,j) = B(iPiv,j)
            B(iPiv,j) = x
          End Do
          T=A(I,I)
        End If
        Do J=1,n
          A(I,J)=A(I,J)/T
          B(I,J)=B(I,J)/T
        End Do
        Do K=1,n
          If (K.Ne.I) Then
            T=A(K,I)
            Do J=1,n
              A(K,J)=A(K,J)-T*A(I,J)
              B(K,J)=B(K,J)-T*B(I,J)
            End Do
          End If
        End Do
      End Do
      DeAllocate ( A  )
      Return
      End ! MatinvPiv

!C***********************************************************************
      Subroutine PrnSig(IOpt,S,xN1,xN2,xN3,S1,S2,S3,P,Q)
      Implicit Double Precision (A-H,O-Z)
      Dimension S(*),xN1(*),xN2(*),xN3(*)

      If (iOpt.Eq.1) Then
        Call Eig_3(0,S,xN1,xN2,xN3,S1,S2,S3,P,Q) ! with Eigenvectors
      Else
        Call Eig_3a(0,S,S1,S2,S3,P,Q) ! no Eigenvectors
      End If
      Return
      End
!C***********************************************************************
      Subroutine Eig_3(iOpt,St,xN1,xN2,xN3,S1,S2,S3,P,Q)
      Implicit Double Precision (A-H,O-Z)
      Dimension St(6),A(3,3),V(3,3), &
               xN1(3),xN2(3),xN3(3)
      !
      ! Get Eigenvalues/Eigenvectors for 3*3 matrix
      ! Wim Bomhof 15/11/'01
      ! PGB : adaption to Principal stress calculation
      !
      ! Applied on principal stresses, directions
      ! Stress vector St(): XX, YY, ZZ, XY, YZ, ZX
      !
      A(1,1) = St(1) ! xx
      A(1,2) = St(4) ! xy = yx
      A(1,3) = St(6) ! zx = xz

      A(2,1) = St(4) ! xy = yx
      A(2,2) = St(2) ! yy
      A(2,3) = St(5) ! zy = yz

      A(3,1) = St(6) ! zx = xz
      A(3,2) = St(5) ! zy = yz
      A(3,3) = St(3) ! zz

      ! Set V to unity matrix
      V(1,1) = 1
      V(2,1) = 0
      V(3,1) = 0

      V(1,2) = 0
      V(2,2) = 1
      V(3,2) = 0

      V(1,3) = 0
      V(2,3) = 0
      V(3,3) = 1


      abs_max_s=0.0
      Do i=1,3
        Do j=1,3
          if (abs(a(i,j)) .Gt. abs_max_s) abs_max_s=abs(a(i,j))
        End Do
      End Do
      Tol = 1d-20 * abs_max_s
      it = 0
      itmax = 50
      Do While ( it.Lt.itMax .And. &
                abs(a(1,2))+abs(a(2,3))+abs(a(1,3)) .Gt. Tol )
        it=it+1
        Do k=1,3
          If (k .Eq. 1) Then
            ip=1
            iq=2
          Else If (k .Eq.2) Then
            ip=2
            iq=3
          Else
            ip=1
            iq=3
          End If
          If (a(ip,iq) .Ne. 0.0) Then
            tau=(a(iq,iq)-a(ip,ip))/(2.0*a(ip,iq))
            If (tau .Ge.0.0) Then
              sign_tau=1.0
            Else
              sign_tau=-1.0
            End If
            t=sign_tau/(abs(tau)+sqrt(1.0+tau*tau))
            c=1.0/sqrt(1.0+t*t)
            s=t*c
            a1p=c*a(1,ip)-s*a(1,iq)
            a2p=c*a(2,ip)-s*a(2,iq)
            a3p=c*a(3,ip)-s*a(3,iq)
            a(1,iq)=s*a(1,ip)+c*a(1,iq)
            a(2,iq)=s*a(2,ip)+c*a(2,iq)
            a(3,iq)=s*a(3,ip)+c*a(3,iq)
            a(1,ip)=a1p
            a(2,ip)=a2p
            a(3,ip)=a3p

            v1p=c*v(1,ip)-s*v(1,iq)
            v2p=c*v(2,ip)-s*v(2,iq)
            v3p=c*v(3,ip)-s*v(3,iq)
            v(1,iq)=s*v(1,ip)+c*v(1,iq)
            v(2,iq)=s*v(2,ip)+c*v(2,iq)
            v(3,iq)=s*v(3,ip)+c*v(3,iq)
            v(1,ip)=v1p
            v(2,ip)=v2p
            v(3,ip)=v3p

            ap1=c*a(ip,1)-s*a(iq,1)
            ap2=c*a(ip,2)-s*a(iq,2)
            ap3=c*a(ip,3)-s*a(iq,3)
            a(iq,1)=s*a(ip,1)+c*a(iq,1)
            a(iq,2)=s*a(ip,2)+c*a(iq,2)
            a(iq,3)=s*a(ip,3)+c*a(iq,3)
            a(ip,1)=ap1
            a(ip,2)=ap2
            a(ip,3)=ap3
          End If ! a(ip,iq)<>0
        End Do ! k
      End Do ! While
      ! principal values on diagonal of a
      S1 = a(1,1)
      S2 = a(2,2)
      S3 = a(3,3)
      ! Derived invariants
      P = (S1+S2+S3)/3
      Q = Sqrt( ( (S1-S2)**2 + (S2-S3)**2 + (S3-S1)**2 ) / 2 )

      ! Sort eigenvalues S1 <= S2 <= S3
      is1 = 1
      is2 = 2
      is3 = 3
      if (s1.Gt.s2) Then
        t   = s2
        s2  = s1
        s1  = t
        it  = is2
        is2 = is1
        is1 = it
      End If
      if (s2.Gt.s3) Then
        t   = s3
        s3  = s2
        s2  = t
        it  = is3
        is3 = is2
        is2 = it
      End If
      if (s1.Gt.s2) Then
        t   = s2
        s2  = s1
        s1  = t
        it  = is2
        is2 = is1
        is1 = it
      End If
      Do i=1,3
        xN1(i) = v(i,is1) ! first  column
        xN2(i) = v(i,is2) ! second column
        xN3(i) = v(i,is3) ! third  column
      End Do
      Return
      End ! Eig_3

      Subroutine Eig_3a(iOpt,St,S1,S2,S3,P,Q) ! xN1,xN2,xN3,
      Implicit Double Precision (A-H,O-Z)
      Dimension St(6),A(3,3)   !  V(3,3),xN1(3),xN2(3),xN3(3)
      !
      ! Get Eigenvalues ( no Eigenvectors) for 3*3 matrix
      ! Wim Bomhof 15/11/'01
      !
      ! Applied on principal stresses, directions
      ! Stress vector XX, YY, ZZ, XY, YZ, ZX
      !
      A(1,1) = St(1) ! xx
      A(1,2) = St(4) ! xy = yx
      A(1,3) = St(6) ! zx = xz

      A(2,1) = St(4) ! xy = yx
      A(2,2) = St(2) ! yy
      A(2,3) = St(5) ! zy = yz

      A(3,1) = St(6) ! zx = xz
      A(3,2) = St(5) ! zy = yz
      A(3,3) = St(3) ! zz

      abs_max_s=0.0
      Do i=1,3
        Do j=1,3
          if (abs(a(i,j)) .Gt. abs_max_s) abs_max_s=abs(a(i,j))
        End Do
      End Do
      Tol = 1d-20 * abs_max_s
      If (iOpt.Eq.1) Tol = 1d-50*abs_max_s
      it=0
      itmax = 50
      Do While ( it.lt.itmax .And.&
                abs(a(1,2))+abs(a(2,3))+abs(a(1,3)) .Gt. Tol )

        it=it+1
        Do k=1,3
          If (k .Eq. 1) Then
            ip=1
            iq=2
          Else If (k .Eq.2) Then
            ip=2
            iq=3
          Else
            ip=1
            iq=3
          End If
          If (a(ip,iq) .Ne. 0.0) Then         ! ongelijk nul ?
            tau=(a(iq,iq)-a(ip,ip))/(2.0*a(ip,iq))
            If (tau .Ge.0.0) Then
              sign_tau=1.0
            Else
              sign_tau=-1.0
            End If
            t=sign_tau/(abs(tau)+sqrt(1.0+tau*tau))
            c=1.0/sqrt(1.0+t*t)
            s=t*c
            a1p=c*a(1,ip)-s*a(1,iq)
            a2p=c*a(2,ip)-s*a(2,iq)
            a3p=c*a(3,ip)-s*a(3,iq)
            a(1,iq)=s*a(1,ip)+c*a(1,iq)
            a(2,iq)=s*a(2,ip)+c*a(2,iq)
            a(3,iq)=s*a(3,ip)+c*a(3,iq)
            a(1,ip)=a1p
            a(2,ip)=a2p
            a(3,ip)=a3p

            ap1=c*a(ip,1)-s*a(iq,1)
            ap2=c*a(ip,2)-s*a(iq,2)
            ap3=c*a(ip,3)-s*a(iq,3)
            a(iq,1)=s*a(ip,1)+c*a(iq,1)
            a(iq,2)=s*a(ip,2)+c*a(iq,2)
            a(iq,3)=s*a(ip,3)+c*a(iq,3)
            a(ip,1)=ap1
            a(ip,2)=ap2
            a(ip,3)=ap3
          End If ! a(ip,iq)<>0
        End Do ! k
      End Do ! While
      ! principal values on diagonal of a
      S1 = a(1,1)
      S2 = a(2,2)
      S3 = a(3,3)
      ! Derived invariants
      P = (S1+S2+S3)/3
      Q = Sqrt( ( (S1-S2)**2 + (S2-S3)**2 + (S3-S1)**2 ) / 2 )

      if (s1.Gt.s2) Then
        t   = s2
        s2  = s1
        s1  = t
      End If
      if (s2.Gt.s3) Then
        t   = s3
        s3  = s2
        s2  = t
      End If
      if (s1.Gt.s2) Then
        t   = s2
        s2  = s1
        s1  = t
      End If
      Return
      End ! Eig_3a

!C
!C***********************************************************************
      Logical Function LEqual(A,B,Eps)
!C***********************************************************************
!C
!C     Returns .TRUE.  when two real values are (almost) equal,
!C             .FALSE. otherwise
!!C
!C I   A,B  : Two real values to be compared
!C I   Eps  : Toleration (Magnitude ~= 1E-5)
!C***********************************************************************
      Implicit Double Precision (A-H,O-Z)
!C***********************************************************************
      LEqual =.True.
      If (A .Eq. B) Return
      If (DAbs(A-B) .LT. 0.5D0*Eps*( DAbs(A) + DAbs(B) + Eps ) )Return
      LEqual =.False.
      Return
      End     ! function LEqual
!C
!C***********************************************************************
      Subroutine CrossProd(xN1,xN2,xN3)
!C***********************************************************************
!!C
!C     Returns cross product of xN1 and xN2
!C
!C I   xN1,xN2 : Two basic vectors
!C O   xN3     : Resulting vector
!C***********************************************************************
      Implicit Double Precision (A-H,O-Z)
      Dimension xN1(*),xN2(*),xN3(*)
!C***********************************************************************

      xN3(1) = xN1(2)*xN2(3) - xN1(3)*xN2(2)
      xN3(2) = xN1(3)*xN2(1) - xN1(1)*xN2(3)
      xN3(3) = xN1(1)*xN2(2) - xN1(2)*xN2(1)

      Return
      End     ! Subroutine CrossProd
!C
!C***********************************************************************
      Double Precision Function ArcSin(X,ie)
!C***********************************************************************
!C
!C     Returns the Arc Sine of X
!C
!C I   X : Input value
!C
!C     Note : In stead of using default routine DASIN we use this one
!C            because ³X³ can be slightly beyond 1 and this will give
!C            a RTE using DASIN(X)
!C
!C***********************************************************************
      Implicit Double Precision (A-H,O-Z)
!C***********************************************************************
      Ie=0
      S = (1-X*X)
!      If (S .Lt. -1E-10) Ie=1
!      If (S .Lt. -1E-10) Write(*,1) X,S
!      If (S .Lt. -1E-10) Write(2,1) X,S
    1 Format(' ArcSin(',1x,1p,e13.5e3,') , S =',1x,1p,e13.5e3)
      If (S.LT.0) S = 0
      S = DSQRT(S)
      ArcSin = DATan2(X,S)
      Return
      End     ! function ArcSin
!C
!C***********************************************************************
      Subroutine CarSig(S1,S2,S3,xN1,xN2,xN3,SNew)
!C***********************************************************************
!C
!C     Returns the Cartesian stresses using the principal stresses S1..S3
!C     and the principal directions
!C
!C I   S1..S3   : Principal stresses
!C I   xN1..xN3 : Principal directions (xNi for Si)
!C
!C***********************************************************************
      Implicit Double Precision (A-H,O-Z)
      Dimension xN1(*),xN2(*),xN3(*),SNew(*)
      Dimension SM(3,3),T(3,3),TT(3,3),STT(3,3)
!C***********************************************************************
!C
!C**** Fill transformation (rotation) matrix
!C
      Do I=1,3
        T(I,1) = xN1(I)
        T(I,2) = xN2(I)
        T(I,3) = xN3(I)
        TT(1,I) = T(I,1)
        TT(2,I) = T(I,2)
        TT(3,I) = T(I,3)
      End Do
!      Call MatTranspose(T,3,TT,3,3,3)

      Call MZeroR(SM,9)
      SM(1,1) = S1
      SM(2,2) = S2
      SM(3,3) = S3
!C
!C**** SMnew = T*SM*TT
!C
      Call MatMat(SM ,3,  TT,3 , 3,3,3 ,STT,3)
      Call MatMat( T ,3, STT,3 , 3,3,3 ,SM ,3)
!     Call MatMatSq(3, SM,  TT, STT )   ! STT = SM*TT
!     Call MatMatSq(3,  T, STT, SM  )   ! SM  =  T*STT
!C
!C**** Extract cartesian stress vector from stress matrix
!C
      Do I=1,3
        SNew(I) = SM(I,I)
      End Do
      SNew(4) = SM(2,1)
      SNew(5) = SM(3,2)
      SNew(6) = SM(3,1)

      Return
      End     ! Subroutine CarSig
!C**********************************************************************
      subroutine setveclen(xn,n,xl)
!C**********************************************************************
      Implicit Double Precision (A-H,O-Z)
      Dimension xN(*)
      x=0
      do i=1,n
        x=x+xn(i)**2
      end do
      if (x.Ne.0) Then
        f=xl/sqrt(x)
        do i=1,3
          xn(i)=xn(i)*f
        end do
      end if
      return
      end ! setveclen

!C**********************************************************************
!C End Of file
!C**********************************************************************
_____________________________________________________________________________________________________
 
MYMODEL1.F90  (SUBROUTINE)
_____________________________________________________________________________________________________
     Subroutine MyModel1 ( IDTask, iMod, IsUndr, &
                           iStep, iTer, iEl, Int, &
                           X, Y, Z, &
                           Time0, dTime, &
                           Props, Sig0, Swp0, StVar0, &
                           dEps, D, BulkW, &
                           Sig, Swp, StVar, ipl, &
                           nStat, &
                           NonSym, iStrsDep, iTimeDep, iTang, &
                           iAbort )
!
! Número de tarea IDTask 1 : Inicializar variables
!                        2 : Calcular tensiones
!                        3 : Calcular la matriz de rigidez
!                        4 : Número de variables de estado
!                        5 : Preguntar propiedades de la matriz
!                        6 : Calcular la matriz de rigidez elástica del material
!
! Argumentos:
!          I/O  Type
!  IDTask   I   I    : Ver arriba
!  iMod     I   I    : Número de modelo (1..10)
!  IsUndr   I   I    : =1 para no drenado, =0 en otro caso
!  iStep    I   I    : Número de paso global
!  iter     I   I    : Número de iteración global
!  iel      I   I    : Número de elemento global
!  Int      I   I    : Número de punto de integración global
!  X        I   R    : Posición X del punto de integración
!  Y        I   R    : Posición Y del punto de integración
!  Z        I   R    : Posición Z del punto de integración
!  Time0    I   R    : Tiempo al inicio del paso de cálculo
!  dTime    I   R    : Incremento de tiempo
!  Props    I   R()  : Lista con los parámetros del modelo
!  Sig0     I   R()  : Tensiones al inicio del paso de cálculo
!  Swp0     I   R    : Exceso d mme presión de poros al inicio del paso de cálculo
!  StVar0   I   R()  : Variable de estado al inicio del paso de cálculo
!  dEps     I   R()  : Incremento de deformación
!  D       I/O  R(,) : Matriz de rigidez del material
!  BulkW   I/O  R    : Módulo volumétrico del agua (sólo caso no drenado)
!  Sig      O   R()  : Tensión resultante
!  Swp      O   R    : Exceso de presión de agua resultante
!  StVar    O   R()  : Resulting values state variables
!  ipl      O   I    : Indicador plástico
!  nStat    O   I    : Número de variables de estado
!  NonSym   O   I    : Matriz D no simétrica
!  iStrsDep O   I    : =1 para matriz D dependiente de las tesniones
!  iTimeDep O   I    : =1 para matriz D dependiente del tiempo
!  iAbort   O   I    : =1 para detenr proceso de cálculo
!
      Implicit Double Precision (A-H, O-Z)

      Dimension Props(*), Sig0(*), StVar0(*), dEps(*), D(6,6), &
                Sig(*), StVar(*)
      
!---  Variables locales
      !Prs_E        : Tensiones principales trial (compresión positiva)
      !Prs          : Tensiones principales al final del paso de cálculo (compresión positiva)
      !Si           : Tensiones principales trial (compresión negativa)
      !xN(i)        : Autocvectores de SigE (xx-yy-zz) = Prs_E (1-2-3)
      
      Dimension dSig(6), xN1(3), xN2(3), xN3(3), xN10(3), xN20(3), xN30(3), &
                Prs_E(3), Prs(3), Prs_0(3)
      
      pi        = 3.1415926535d0
      Rad       = pi/180d0
      io        = 0
      ttol      = 1d-9
      nStatV    = 2
      ideb      = 1 !para debugging

!Contenido de Props()
!                I/O  Type      
! 1:   G          I    R     Módulo de corte
! 2:   nu         I    R     Coeficiente de Poisson 
! 3:   E50_ref    I    R     Rigidez triaxial secante (superficie cono)
! 4:   Eoed_ref   I    R     Rigidez edométrica tangente (superficie cap)
! 5:   Eur_ref    I    R     Rigidez en descarga/recarga
! 6:   xm         I    R     Potencia de la dependencia tensional
! 7:   c          I    R     Cohesión (efectva)
! 8:   phi        I    R     Ángulo de fricción (efectiva)
! 9:   psi        I    R     Ángulo de dilatancia
! 10:  xnu_ur     I    R     Coeficiente de poisson en descarga/recarga
! 11:  p_ref      I    R     Presión de refencia para la rigidez
! 12:  xK0_nc     I    R     K0 (normalmente consolidado)
! 13:  R_f        I    R     Ratio de rotura
! 14:  Sigma_trac I    R     Máxima tensión de tracción (aplica sobre Prs, tracción +, compresión -)
      
! Prs_E           I    R(3) : Tensiones principales del predictor elástico S1 >= S2 >= S3
! Prs             O    R(3) : Tensiones principales resultantes 
! ipl             O    I    : Indicador de punto plástico
!                             1 : Superfície límite de Coulomb (MC_plast)
!                             2 : Superficie de tracción (vertice)
!                             3 : Punto plástico en CAP (CAP_plast)
!                             4 : Punto plástico en HS_CAP (HS_CAP_plast)
!                             5 : Punto plástico en HS (HS_plast)

!Contenido de StVar0
!                I/O  Type            
! 1:  gamma_p    I/O   R     Variable de endurecimiento HS
! 2:  p_p        I/O   R     Variable de endurecimiento CAP
      
      If (IDTask .EQ. 1) Then !Inicializo variables

        !Parámetros y variables de estado

        G          = Props(1) 
        xNu        = Props(2)     
        E50_ref    = Props(3) 
        Eoed_ref   = Props(4)
        Eur_ref    = Props(5)
        xm         = Props(6)
        c          = Props(7)
        phi        = Props(8)*Rad
        psi        = Props(9)*Rad
        xNu_ur     = Props(10)
        p_ref      = Props(11)
        xK0_nc     = Props(12)
        R_f        = Props(13)
        sigma_trac = Props(14)

        If (iStep .EQ. 0) Then   

        !Cálculo variables endurecimiento en primer passo de cálculo
        iOpt        = 0
        Call PrnSig(iOpt, Sig0, xN10, xN20, xN30, S10, S20, S30, P, Q)
        Prs_0(1)    = S10*(-1d0)
        Prs_0(2)    = S20*(-1d0)
        Prs_0(3)    = S30*(-1d0)
        
        !Superficie HS (StVar0(1))
        qf1         = max(c*cos(phi)/sin(phi)+Prs_0(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))
        qf2         = (2d0*sin(phi)/(1d0-sin(phi)))
        qf          = qf1*qf2
        qa          = qf/R_f
        E_fac1      = max(c*cos(phi)+Prs_0(3)*sin(phi), max(c*cos(phi)-sigma_trac*sin(phi), 0d0))
        E_fac2      = (c*cos(phi)+p_ref*sin(phi))
        E50         = E50_ref*(E_fac1/E_fac2)**xm
        Ei          = 2d0*E50/(2d0-R_f)  
        Eur         = Eur_ref*(E_fac1/E_fac2)**xm           
        If (abs(Prs_E(1)-Prs_E(3)) .LT. qf) Then
        f13_hs_fac1 = ((Prs_0(1)-Prs_0(3))/(qa-(Prs_0(1)-Prs_0(3))))
        f13_hs_fac2 = 2d0*(Prs_0(1)-Prs_0(3))/Eur
        gamma_p_0   = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2            
        Else
        f13_hs_fac1 = qf/(qa-qf)
        f13_hs_fac2 = 2d0*qf/Eur
        gamma_p_0   = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2            
        End If
        
        !Superficie CAP (StVar0(2))
        xlam_ast    = p_ref/Eoed_ref
        xkap_ast    = 3d0*p_ref*(1-2*xNu_ur)/Eur_ref
        alfa1       = ((1d0-xK0_nc)**2d0)/((1d0+2d0*xK0_nc)**2d0)
        alfa21      = (1d0-xK0_nc)*(1d0-2d0*xNu_ur)*(xlam_ast/xkap_ast-1d0)
        alfa221     = (1d0+2d0*xK0_nc)*(1d0-2d0*xNu_ur)*(xlam_ast/xkap_ast)
        alfa222     = (1d0-xK0_nc)*(1d0+xNu_ur)
        alfa        = 3d0*(alfa1+(alfa21/(alfa221-alfa222)))**0.5d0              
        delta       = (3d0+sin(phi))/(3d0-sin(phi))
        q_m_0       = Prs_0(1)+(delta-1d0)*Prs_0(2)-delta*Prs_0(3)
        p_0         = (Prs_0(1)+Prs_0(2)+Prs_0(3))/3d0
        p_p_0       = ((q_m_0**2d0)/(alfa**2d0)+p_0**2d0)**0.5d0
        
        !Cálculo variables en primer paso de cálculo
        StVar0(1)   = max(max(StVar0(1), gamma_p_0), 1d-10)
        StVar0(2)   = min(StVar0(2), min((-1d0)*p_p_0,(-1d0)*c*cos(phi)/sin(phi)))
        End If
        
        !Corrección valor xNu 
        If (xNu .GT. 0.5d0-1d-6 .OR. xNu .LT. -1d-6) Then
          Call OK_MessageBox('Invalid value for "nu" (param 2)')
          iAbort    = 1
          Return
        End If
        If (xNu .LT. 0d0) Then
          xNu       = 0.0d0
          Props(2)  = xNu
        End If
        If (xNu .GT. 0.495d0) Then
          xNu       = 0.495d0
          Props(2)  = xNu
        End If
        Call CopyRVec(StVar0, StVar, nStatV) !Copio StVar0 en StVar  
      
      End If  ! IDTask = 1

    
      If (IDTask .EQ. 2) Then !Cálculo de tensiones

      Write(1,*) 'IDTask:', IDTask, 'iEl:', iEl, 'Int:', Int, 'iStep:', iStep, 'iTer:', iTer               
      Call WriIvl(ideb, 'IsUndr', IsUndr)

        !Parámetros y variables de estado
        Call CopyRVec(StVar0, StVar, nStatV) !Copio StVar0 en StVar
        G         = Props(1) 
        xNu       = Props(2)     
        E50_ref   = Props(3) 
        Eoed_ref  = Props(4)
        Eur_ref   = Props(5)
        xm        = Props(6)
        c         = Props(7)
        phi       = Props(8)*Rad
        psi       = Props(9)*Rad
        xNu_ur    = Props(10)
        p_ref     = Props(11)
        xK0_nc    = Props(12)
        R_f       = Props(13)
        sigma_trac = Props(14)   
        
        gamma_p_E = StVar0(1)
        p_p_E     = (-1d0)*StVar0(2)
        
        Call WriVal(ideb, 'gamma_p_E', gamma_p_E)
        Call WriVal(ideb, 'p_p_E', p_p_E)
        
        !Construcción matriz D6
        iOpt        = 0
        Call PrnSig(iOpt, Sig0, xN10, xN20, xN30, S10, S20, S30, P0, Q0)
        Prs_0(1)    = S10*(-1d0)
        Prs_0(2)    = S20*(-1d0)
        Prs_0(3)    = S30*(-1d0)        
        Call WriVec(ideb, 'Prs_0', Prs_0, 3)
        Call Matriz6(Props, Prs_0, D, IsUndr)
        Call WriMat(ideb, 'D6' , D , 6, 6, 6)
                
        !Paso de cálculo elástico
        Call MatVec(D, 6, dEps, 6, dSig)
        Call AddVec(Sig0, dSig, 1d0, 1d0, 6, Sig )      

        Call WriVec(ideb,'dEps', dEps, 6)              
        Call WriVec(ideb,'Sig0', Sig0,6)                     
        Call WriVec(ideb,'Sig', Sig,6)    
        
        iOpt        = 1
        Call PrnSig(iOpt, Sig, xN1, xN2, xN3, S1, S2, S3, P, Q)
        !Sig       :  compresión negativa (S1<S2<S3)
        !Prs, Prs_E:  compresión positiva (S1>S2>S3)
        Prs_E(1)    = S1*(-1d0)
        Prs_E(2)    = S2*(-1d0)
        Prs_E(3)    = S3*(-1d0)
        p_E         = (Prs_E(1)+Prs_E(2)+Prs_E(3))/3d0
        
        Call WriVec(ideb, 'Prs_E', Prs_E, 3)
                
        !Cálculo superficies de fluencia f13_hs, f13_mc y f_cap en Prs_E

        !Cálculo f13_hs en Prs_E 
        qf1         = max(c*cos(phi)/sin(phi)+Prs_E(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))
        qf2         = (2d0*sin(phi)/(1d0-sin(phi)))
        qf_E        = qf1*qf2
        qa          = qf_E/R_f
        E_fac1      = max(c*cos(phi)+Prs_E(3)*sin(phi), max(c*cos(phi)-sigma_trac*sin(phi), 0d0))
        E_fac2      = (c*cos(phi)+p_ref*sin(phi))
        E50         = E50_ref*(E_fac1/E_fac2)**xm
        Ei          = 2d0*E50/(2d0-R_f)  
        Eur         = Eur_ref*(E_fac1/E_fac2)**xm   

        If (abs(Prs_E(1)-Prs_E(3)) .LT. qf_E) Then         
        f13_hs_fac1 = ((Prs_E(1)-Prs_E(3))/(qa-(Prs_E(1)-Prs_E(3))))
        f13_hs_fac2 = 2d0*(Prs_E(1)-Prs_E(3))/Eur
        f13_hs      = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p_E        
        Else
        f13_hs      = 1d0
        End If
        
        Call WriVal(ideb, 'Prs_E(1)-Prs_E(3)', Prs_E(1)-Prs_E(3))                      
        Call WriVal(ideb, 'qa', qa)  
        Call WriVal(ideb, 'qa-(Prs_E(1)-Prs_E(3))', qa-(Prs_E(1)-Prs_E(3)))                  
        Call WriVal(ideb, 'f13_hs', f13_hs)              
        
        !Cálculo f13_mc en PrS_E
        f13_mc_fac1 = (1d0/2d0)*(Prs_E(1)-Prs_E(3))
        f13_mc_fac2 = (1d0/2d0)*(Prs_E(1)+Prs_E(3))*sin(phi)
        f13_mc      = f13_mc_fac1-f13_mc_fac2-c*cos(phi)

        Call WriVal(ideb, 'f13_mc', f13_mc)        
        
        !Cálculo f_cap en Prs_E
        xlam_ast    = p_ref/Eoed_ref
        xkap_ast    = 3d0*p_ref*(1-2*xNu_ur)/Eur_ref
        alfa1       = ((1d0-xK0_nc)**2d0)/((1d0+2d0*xK0_nc)**2d0)
        alfa21      = (1d0-xK0_nc)*(1d0-2d0*xNu_ur)*(xlam_ast/xkap_ast-1d0)
        alfa221     = (1d0+2d0*xK0_nc)*(1d0-2d0*xNu_ur)*(xlam_ast/xkap_ast)
        alfa222     = (1d0-xK0_nc)*(1d0+xNu_ur)
        alfa        = 3d0*(alfa1+(alfa21/(alfa221-alfa222)))**0.5d0              
        delta       = (3d0+sin(phi))/(3d0-sin(phi))
        q_m         = Prs_E(1)+(delta-1d0)*Prs_E(2)-delta*Prs_E(3)
        p           = (Prs_E(1)+Prs_E(2)+Prs_E(3))/3d0      
        f_cap       = (q_m**2d0)/(alfa**2d0)+(p**2d0)-(p_p_E**2d0)            
        f_cap_E     = f_cap
         
         Call WriVal(ideb, 'f_cap', f_cap)                
 !       Call WriVal(ideb, 'xK0_nc', xK0_nc)   
 !       Call WriVal(ideb, 'alfa', alfa)        
 !       Call WriVal(ideb, 'xH', xH)          
        
        !Cálculo f_trac en Prs_E        
        f13_trac = sigma_trac-Prs_E(3)
        
        Call WriVal(ideb, 'f13_trac', f13_trac)        
        
        !Cálculo tolerancias 
        f13_hs_TOL    = ttol*abs(f13_hs)+ttol
        f13_mc_TOL    = ttol*abs(f13_mc)+ttol
        f_cap_TOL     = ttol*abs(f_cap)+ttol
        f13_trac_TOL  = ttol*abs(f13_trac)+ttol
            
        !Call WriVal(ideb, 'f13_hs_TOL', f13_hs_TOL)       
        !Call WriVal(ideb, 'f13_mc_TOL', f13_mc_TOL)       
        !Call WriVal(ideb, 'f_cap_TOL', f_cap_TOL)            
        !Call WriVal(ideb, 'f13_trac_TOL', f13_trac_TOL)                            
        
        !Comprobación si el paso de scálculo es elástico
        If (f13_hs .LT. f13_hs_TOL .AND. &
            f13_mc .LT. f13_mc_TOL .AND. &            
            f13_trac .LT. f13_trac_TOL .AND. &                        
            f_cap .LT. f_cap_TOL) Then        
            write(1,*) "PASO DE CÁLCULO ELÁSTICO, DENTRO DE f13_hs, f13_mc, f_cap y f13_trac"
            Prs(1)     = Prs_E(1)
            Prs(2)     = Prs_E(2)            
            Prs(3)     = Prs_E(3)    
            ipl        = 0
            GOTO 800    
        End If         
           
        !Comprobación si Prs_E está fuera de f13_hs
        If (f13_hs .GE. f13_hs_TOL) Then 
        If (abs(Prs_E(1)-Prs_E(3)) .LT. qf_E) Then 
        write(1,*) "ESTOY FUERA DE f13_hs"
        !Cálculo Prs en superficies de HS y endurecimiento gamma_p (StVar)
          iArea=1 !Calculo en plano f13_hs
          Call HS_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr) 
          !Comprobación cambios en el orden de las tensiones principales
          Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN HS', Prs, 3)          
          If (Prs(2) .LT. Prs(3)) Then 
          write(1,*) "ESTOY FUERA DE f13_hs, AREA=2"              
          iArea = 2 !Tr. Compresion (estoy en eje S2=S3)  
          End If
          If (Prs(1) .LT. Prs(2)) Then 
          write(1,*) "ESTOY FUERA DE f13_hs, AREA=3"                            
          iArea = 3 !Tr. Extension (estoy en eje S1=S2)              
          End If
          If (iArea .NE. 1) Then !Returning Map en eje 
          Call HS_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
          Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN HS, AREA 2,3', Prs, 3)          
          End If
        Else 
        Prs(1)      = Prs_E(1)
        Prs(2)      = Prs_E(2)
        Prs(3)      = Prs_E(3)
        End If
                          
          !Comprobación si nuevas tensiones Prs no superan MC
          f13_mc_fac1 = (1d0/2d0)*(Prs(1)-Prs(3))
          f13_mc_fac2 = (1d0/2d0)*(Prs(1)+Prs(3))*sin(phi)
          f13_mc      = f13_mc_fac1-f13_mc_fac2-c*cos(phi)

          Call WriVal(ideb, 'f13_mc', f13_mc) 
          
            If (f13_mc .GE. f13_mc_TOL) Then 
            !No sirve Prs anterior
            !Desendurecimiento HS adoptando StVar0 y se vuelve a aplicar a Prs_E de partida
            !Cálculo Prs en superficies MC (no tienen endurecimiento)
            iArea = 1
            write(1,*) "ESTOY FUERA DE f13_mc"
            Call MC_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)                      
            !Comprobación cambios en el orden de las tensiones principales
            Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN MC', Prs, 3)                      
            If (Prs(2) .LT. Prs(3)) Then 
            write(1,*) "ESTOY FUERA DE f13_mc, AREA=2" 
            iArea = 2 !Tr. Compresion (estoy en eje S2=S3)  
            End If
            If (Prs(1) .LT. Prs(2)) Then 
            write(1,*) "ESTOY FUERA DE f13_mc, AREA=3"                            
            iArea = 3 !Tr. Extension (estoy en eje S1=S2)              
            End If
            If (iArea .NE. 1) Then !Returning Map en eje 
            Call MC_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)                     
            Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN MC, AREA 2,3', Prs, 3)                      
            End If
            End If
            
          !Comprobación si nuevas tensiones Prs no superan TRAC
            f13_trac = sigma_trac-Prs(3)
            
            Call WriVal(ideb, 'f13_trac', f13_trac)        
            
            If (f13_trac .GE. f13_trac_TOL) Then 
            !No sirve Prs anterior
            !Desendurecimiento HS adoptando StVar0 y se vuelve a aplicar a Prs_E de partida
            !Cálculo Prs en superficie TRAC (no tienen endurecimiento)
            iArea = 1
            write(1,*) "ESTOY FUERA DE f13_trac"              
            Call TRAC_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)         
            !Comprobación cambios en el orden de las tensiones principales
            Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN TRAC', Prs, 3)                        
            If (Prs(2) .LT. Prs(3)) Then 
            write(1,*) "ESTOY FUERA DE f13_trac, AREA=2"              
            iArea = 2 !Tr. Compresion (eje S2=S3)  
            End If
            If (Prs(1) .LT. Prs(2)) Then 
            write(1,*) "ESTOY FUERA DE f12, AREA=3"              
            iArea = 3 !Tr. Extension (eje S1=S2)              
            End If
            If (iArea .NE. 1) Then !Returning Map en eje 
            Call TRAC_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
            Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN TRAC, AREA 2,3', Prs, 3)                      
            End If
         
            !Comprobación si nuevas tensiones Prs no superan HS            
            qf1         = max(c*cos(phi)/sin(phi)+Prs(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))
            qf2         = (2d0*sin(phi)/(1d0-sin(phi)))
            qf          = qf1*qf2
            qa          = qf/R_f
            E_fac1      = max(c*cos(phi)+Prs(3)*sin(phi), max(c*cos(phi)-sigma_trac*sin(phi), 0d0))
            E_fac2      = (c*cos(phi)+p_ref*sin(phi))
            Eur         = Eur_ref*(E_fac1/E_fac2)**xm
            E50         = E50_ref*(E_fac1/E_fac2)**xm
            Ei          = 2d0*E50/(2d0-R_f)     
            
            If (abs(Prs(1)-Prs(3)) .LT. qf) Then                 
            f13_hs_fac1 = ((Prs(1)-Prs(3))/(qa-(Prs(1)-Prs(3))))
            f13_hs_fac2 = 2d0*(Prs(1)-Prs(3))/Eur
            f13_hs      = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p_E             
            Else
            f13_hs      = 1d0
            End  If
             
              If (f13_hs .GE. f13_hs_TOL) Then            
              If (abs(Prs_E(1)-Prs_E(3)) .LT. qf_E) Then 
              !No sirve Prs anterior
              !Desendurecimiento HS adoptando StVar0 y se vuelve a aplicar a Prs_E de partida
              !Cálculo Prs en superficies HS y TRAC (no tienen endurecimiento)
              iArea = 1
              write(1,*) "ESTOY FUERA DE f13_hs y f13_trac"
              Call HS_TRAC_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)                     
              !Comprobación cambios en el orden de las tensiones principales
              Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN HS Y TRAC', Prs, 3)                      
              If (Prs(2) .LT. Prs(3)) Then 
              write(1,*) "ESTOY FUERA DE f13_hs y f13_trac, AREA=2" 
              iArea = 2 !Tr. Compresion (estoy en eje S2=S3)  
              End If
              If (Prs(1) .LT. Prs(2)) Then 
              write(1,*) "ESTOY FUERA DE f13_hs y f13_trac, AREA=3"                            
              iArea = 3 !Tr. Extension (estoy en eje S1=S2)              
              End If
              If (iArea .NE. 1) Then !Returning Map en eje 
              Call HS_TRAC_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
              Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN HS Y TRAC, AREA 2,3', Prs, 3)                                
              End If
              End If
              
                !Comprobación si nuevas tensiones Prs no superan MC
                f13_mc_fac1 = (1d0/2d0)*(Prs(1)-Prs(3))
                f13_mc_fac2 = (1d0/2d0)*(Prs(1)+Prs(3))*sin(phi)
                f13_mc      = f13_mc_fac1-f13_mc_fac2-c*cos(phi)
                
                Call WriVal(ideb, 'f13_mc', f13_mc) 
            
                If (f13_mc .GE. f13_mc_TOL) Then 
                !No sirve Prs anterior
                !Desendurecimiento HS adoptando StVar0 y se vuelve a aplicar a Prs_E de partida
                !Cálculo Prs en superficies MC y TRAC (no tienen endurecimiento)
                iArea = 1
                write(1,*) "ESTOY FUERA DE f13_mc y f13_trac"
                Call MC_TRAC_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)         
                !Comprobación cambios en el orden de las tensiones principales
                Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN MC y TRAC', Prs, 3)                      
                If (Prs(2) .LT. Prs(3)) Then 
                iArea = 2 !Tr. Compresion (estoy en eje S2=S3)  
                End If
                If (Prs(1) .LT. Prs(2)) Then 
                write(1,*) "ESTOY FUERA DE f13_mc y f13_trac, AREA=3"                            
                iArea = 3 !Tr. Extension (estoy en eje S1=S2)              
                End If
                If (iArea .NE. 1) Then !Returning Map en eje 
                Call MC_TRAC_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)         
                Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN MC y TRAC, AREA 2,3', Prs, 3)                      
                End If
                End If            
           End If
           End If 
          
          !Compruebo que las nuevas tensiones Prs (de HS o MC) no superan f_cap (alfa calculado antes)  
          q_m         = Prs(1)+(delta-1d0)*Prs(2)-delta*Prs(3)
          p           = (Prs(1)+Prs(2)+Prs(3))/3d0          
          f_cap       = (q_m**2d0)/(alfa**2d0)+p**2d0-p_p_E**2d0            

          Call WriVal(ideb, 'f_cap', f_cap) 
          
          If (f_cap .GE. f_cap_TOL) Then 
          !No sirve Prs anterior
          !Desendurecimiento HS adoptando StVar0 y se vuelve a aplicar a Prs_E de partida
          !Calculo Prs en superficie CAP y endurezco p_p
            write(1,*) "ESTOY FUERA DE f_cap"  
            If (f_cap_E .GT. f_cap_TOL) Then
            Call CAP_plast(Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)                     
            End If                
            !Compruebo que las nuevas tensiones Prs no superan f13_hs                        
            qf1         = max(c*cos(phi)/sin(phi)+Prs(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))
            qf2         = (2d0*sin(phi)/(1d0-sin(phi)))
            qf          = qf1*qf2
            qa          = qf/R_f
            E_fac1      = max(c*cos(phi)+Prs(3)*sin(phi), max(c*cos(phi)-sigma_trac*sin(phi), 0d0))
            E_fac2      = (c*cos(phi)+p_ref*sin(phi))
            Eur         = Eur_ref*(E_fac1/E_fac2)**xm
            E50         = E50_ref*(E_fac1/E_fac2)**xm
            Ei          = 2d0*E50/(2d0-R_f)     
        
            If (abs(Prs(1)-Prs(3)) .LT. qf) Then                 
            f13_hs_fac1 = ((Prs(1)-Prs(3))/(qa-(Prs(1)-Prs(3))))
            f13_hs_fac2 = 2d0*(Prs(1)-Prs(3))/Eur
            f13_hs      = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p_E
            Else 
            f13_hs      = 1d0
            End If

            Call WriVal(ideb, 'f13_hs', f13_hs) 
            
              If (f13_hs .GE. f13_hs_TOL) Then        
              If (abs(Prs_E(1)-Prs_E(3)) .LT. qf_E) Then 
              !Mixed Return HS-CAP
              !No sirve Prs anterior
              !Desendurecimiento CAP adoptando StVar0 y se vuelve a aplicar a Prs_E de partida
              !Calculo Prs en superficie HS y CAP y endurecimiento gamma_p y p_p
              iArea = 1
              write(1,*) "ESTOY FUERA DE f13_hs y f_cap"              
              Call HS_CAP_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
              !Comprobación cambios en el orden de las tensiones principales
              Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN HS y CAP', Prs, 3)                        
              If (Prs(2) .LT. Prs(3)) Then 
              write(1,*) "ESTOY FUERA DE f13_hs y f_cap, AREA=2"              
              iArea = 2 !Tr. Compresion (eje S2=S3)  
              End If
              If (Prs(1) .LT. Prs(2)) Then 
              write(1,*) "ESTOY FUERA DE f13_hs y f_cap, AREA=3"              
              iArea = 3 !Tr. Extension (eje S1=S2)              
              End If
              If (iArea .NE. 1) Then !Returning Map en eje 
              Call HS_CAP_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
              Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN HS y CAP, AREA 2,3', Prs, 3)                        
              End If
              End If
            
              !Compruebo que las nuevas tensiones Prs no superan MC
              f13_mc_fac1 = (1d0/2d0)*(Prs(1)-Prs(3))
              f13_mc_fac2 = (1d0/2d0)*(Prs(1)+Prs(3))*sin(phi)
              f13_mc      = f13_mc_fac1-f13_mc_fac2-c*cos(phi)

              Call WriVal(ideb, 'f13_mc', f13_mc)               
                          
              If (f13_mc .GE. f13_mc_TOL) Then 
              !Returning Map mixto MC-CAP
              !No sirve Prs anterior
              !Desendurecimiento HS y CAP adoptando StVar0 y se vuelve a aplicar a Prs_E de partida
              !Calculo Prs en superficies CAP y MC y endurecimiento gamma_p 
              iArea = 1
              write(1,*) "ESTOY FUERA DE f13_mc y f_cap"                            
              Call MC_CAP_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)              
              !Comprobación cambios en el orden de las tensiones principales
              Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN MC y CAP', Prs, 3)                        
              If (Prs(2) .LT. Prs(3)) Then 
              write(1,*) "ESTOY FUERA DE f13_mc y f_cap, AREA=2"              
              iArea = 2 !Tr. Compresion (eje S2=S3)  
              End If
              If (Prs(1) .LT. Prs(2)) Then 
              write(1,*) "ESTOY FUERA DE f13_mc y f_cap, AREA=3"                                
              iArea = 3 !Tr. Extension (eje S1=S2)              
              End If
              If (iArea .NE. 1) Then !Returning Map en eje 
              Call MC_CAP_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)         
              Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN MC y CAP, AREA 2,3', Prs, 3)                        
              End If
              End If              
            End If            
          End If 
          
        Else 
                      
        !Compruebo si Prs_E está fuera de f_cap (alfa calculado antes)
        q_m       = Prs_E(1)+(delta-1)*Prs_E(2)-delta*Prs_E(3)
        p         = (Prs_E(1)+Prs_E(2)+Prs_E(3))/3d0            
        f_cap     = (q_m**2d0)/(alfa**2d0)+p**2d0-p_p_E**2d0            
        
          If (f_cap .GE. f_cap_TOL) Then         
          !Calculo Prs en superficie CAP y endurezco p_p
            write(1,*) "ESTOY FUERA DE f_cap"                
            Call CAP_plast(Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)         
            !Compruebo que las nuevas tensiones Prs no superan f13_hs            
            qf1         = max(c*cos(phi)/sin(phi)+Prs(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))
            qf2         = (2d0*sin(phi)/(1d0-sin(phi)))
            qf          = qf1*qf2
            qa          = qf/R_f
            E_fac1      = max(c*cos(phi)+Prs(3)*sin(phi), max(c*cos(phi)-sigma_trac*sin(phi), 0d0))
            E_fac2      = (c*cos(phi)+p_ref*sin(phi))
            Eur         = Eur_ref*(E_fac1/E_fac2)**xm
            E50         = E50_ref*(E_fac1/E_fac2)**xm
            Ei          = 2d0*E50/(2d0-R_f)     

            If (abs(Prs(1)-Prs(3)) .LT. qf) Then             
            f13_hs_fac1 = ((Prs(1)-Prs(3))/(qa-(Prs(1)-Prs(3))))
            f13_hs_fac2 = 2d0*(Prs(1)-Prs(3))/Eur
            f13_hs      = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p_E
            Else
            f13_hs      = 1d0
            End If
            
            Call WriVal(ideb, 'f13_hs', f13_hs) 
            
              If (f13_hs .GE. f13_hs_TOL) Then
              If (abs(Prs_E(1)-Prs_E(3)) .LT. qf_E) Then 
              !Returning Map mixto HS-CAP
              !No sirve Prs anterior
              !Desendurecimiento CAP adoptando StVar0 y se vuelve a aplicar a Prs_E de partida
              !Calculo Prs en superficie HS y CAP y endurecimiento gamma_p y p_p
              iArea = 1
              write(1,*) "ESTOY FUERA DE f13_hs y f_cap"              
              Call HS_CAP_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
              !Comprobación cambios en el orden de las tensiones principales
              Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN HS', Prs, 3)                        
              If (Prs(2) .LT. Prs(3)) Then 
              write(1,*) "ESTOY FUERA DE f13_hs y f_cap, AREA=2"              
              iArea = 2 !Tr. Compresion (eje S2=S3)  
              End If
              If (Prs(1) .LT. Prs(2)) Then 
              write(1,*) "ESTOY FUERA DE f13_hs y f_cap, AREA=3"              
              iArea = 3 !Tr. Extension (eje S1=S2)              
              End If
              If (iArea .NE. 1) Then !Returning Map en eje 
              Call HS_CAP_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
              Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN HS, AREA 2,3', Prs, 3)                        
              End If
              End If
              End If
              
              !Compruebo que las nuevas tensiones Prs no superan MC
              f13_mc_fac1 = (1d0/2d0)*(Prs(1)-Prs(3))
              f13_mc_fac2 = (1d0/2d0)*(Prs(1)+Prs(3))*sin(phi)
              f13_mc      = f13_mc_fac1-f13_mc_fac2-c*cos(phi)

              Call WriVal(ideb, 'f13_mc', f13_mc)               
                          
              If (f13_mc .GE. f13_mc_TOL) Then 
              !Returning Map mixto MC-CAP
              !No sirve Prs anterior 
              !Desendurecimiento HS y CAP adoptando StVar0 y se vuelve a aplicar a Prs_E de partida
              !Calculo Prs en superficies CAP y MC y endurecimiento gamma_p 
              iArea = 1
              write(1,*) "ESTOY FUERA DE f13_mc y f_cap"              
              Call MC_CAP_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)       
              !Comprobación cambios en el orden de las tensiones principales
              Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN MC', Prs, 3)                        
              If (Prs(2) .LT. Prs(3)) Then 
              write(1,*) "ESTOY FUERA DE f13_mc y f_cap, AREA=2"              
              iArea = 2 !Tr. Compresion (eje S2=S3)  
              End If
              If (Prs(1) .LT. Prs(2)) Then 
              write(1,*) "ESTOY FUERA DE f13_mc y f_cap, AREA=3"              
              iArea = 3 !Tr. Extension (eje S1=S2)              
              End If
              If (iArea .NE. 1) Then !Returning Map en eje 
              Call MC_CAP_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)            
              Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN MC, AREA 2,3', Prs, 3)                        
              End If
            End If            
          End If
        End If
    
        !Cálculo f13_trac, f13_hs, f_cap y f13_mc con Prs_E
        qf1         = max(c*cos(phi)/sin(phi)+Prs_E(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))
        qf2         = (2d0*sin(phi)/(1d0-sin(phi)))
        qf_E        = qf1*qf2
        qa          = qf_E/R_f
        E_fac1      = max(c*cos(phi)+Prs_E(3)*sin(phi), max(c*cos(phi)-sigma_trac*sin(phi), 0d0))
        E_fac2      = (c*cos(phi)+p_ref*sin(phi))
        Eur         = Eur_ref*(E_fac1/E_fac2)**xm
        E50         = E50_ref*(E_fac1/E_fac2)**xm
        Ei          = 2d0*E50/(2d0-R_f)     
            
        If (abs(Prs_E(1)-Prs_E(3)) .LT. qf_E) Then         
        f13_hs_fac1 = ((Prs_E(1)-Prs_E(3))/(qa-(Prs_E(1)-Prs_E(3))))
        f13_hs_fac2 = 2d0*(Prs_E(1)-Prs_E(3))/Eur
        f13_hs      = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p_E
        Else 
        f13_hs      = 1d0
        End If

        f13_trac = sigma_trac-Prs_E(3)        
        
        If (f13_trac .GE. f13_trac_TOL .AND. &
            f13_hs .LT. f13_hs_TOL) Then 
        !Calculo Prs en superficies TRAC 
        iArea = 1
        write(1,*) "ESTOY FUERA DE f13_trac"              
        Call TRAC_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)         
        !Comprobación cambios en el orden de las tensiones principales
        Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN TRAC', Prs, 3)                        
        If (Prs(2) .LT. Prs(3)) Then 
        write(1,*) "ESTOY FUERA DE f13_trac, AREA=2"              
        iArea = 2 !Tr. Compresion (eje S2=S3)  
        End If
        If (Prs(1) .LT. Prs(2)) Then 
        write(1,*) "ESTOY FUERA DE f12, AREA=3"              
        iArea = 3 !Tr. Extension (eje S1=S2)              
        End If
        If (iArea .NE. 1) Then !Returning Map en eje 
        Call TRAC_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
        Call WriVec(ideb, 'Prs TRAS RETURNING MAP EN TRAC, AREA 2,3', Prs, 3)                        
        End If
        End If              
            
        p = (Prs(1)+Prs(2)+Prs(3))/3d0
        If (p .LT. (-1d0)*sigma_trac) Then  !Si p<-sigma_trac quiere decir que se ha hecho
                                            !returning map en cono complementario
            Prs(1) = (-1d0)*sigma_trac
            Prs(2) = (-1d0)*sigma_trac
            Prs(3) = (-1d0)*sigma_trac
            ipl=2
        End If
            
800     Continue        

        Write(1,*) "Las tensiones ppales al final de este paso de cálculo son Prs="            
        Call WriVec(1, 'Prs', Prs, 3)            

        Call WriVal(ideb, 'gamma_p_E', StVar0(1))
        Call WriVal(ideb, 'gamma_p', StVar(1))
        Call WriVal(ideb, 'p_p_E', (-1d0)*StVar0(2))        
        Call WriVal(ideb, 'p_p', (-1d0)*StVar(2))        
        
        !Se vuelve a las tensiones en coordenadas cartesianas        
        !Prs: tracción negativa, compresión positiva
        S1 = (-1d0)*Prs(1)
        S2 = (-1d0)*Prs(2)
        S3 = (-1d0)*Prs(3)
        Call CarSig(S1, S2, S3, xN1, xN2, xN3, Sig)
        !Sig: tracción positiva, compresión negativa

 !       Write(1,*) "Las tensiones al final de este paso de cálculo son Sig"              
 !       Call WriVec(ideb, 'Sig', Sig, 6)        
        
        E_fac1      = max(c*cos(phi)+Prs_0(3)*sin(phi), max(c*cos(phi)-sigma_trac*sin(phi), 0d0))
        E_fac2      = (c*cos(phi)+p_ref*sin(phi))
        Eur         = Eur_ref*(E_fac1/E_fac2)**xm            
        G           = Eur/(2d0*(1d0+xNu_ur))        
        xNu_U       = 0.495d0
        Fac         = (1d0+xNu_U)/(1d0-2d0*xNu_U)-(1d0+xNu_ur)/(1d0-2d0*xNu_ur)
        
        If (IsUndr .EQ. 1) Then
        BulkW       = 2d0*G/3d0*Fac                
        Else
        BulkW       = 0d0
        End If
        
        !Caso no drenado
        If (IsUndr .EQ. 1) Then
          dEpsV     = dEps(1)+dEps(2)+dEps(3)
          dSwp      = BulkW*dEpsV
          Swp       = Swp0+dSwp
        Else
          Swp       = Swp0
        End If
        
      End If  ! IDTask = 2

      
      If (IDTask .EQ. 3 .OR. &
         IDTask .EQ. 6 ) Then ! Cálculo de la matriz D                
        !Construcción matriz D
        iOpt        = 1
        Call PrnSig(iOpt, Sig0, xN1, xN2, xN3, S10, S20, S30, P0, Q0)
        Prs_0(1)    = S10*(-1d0)
        Prs_0(2)    = S20*(-1d0)
        Prs_0(3)    = S30*(-1d0)        
        Call Matriz6(Props, Prs_0, D, IsUndr)    
        BulkW = 0
        If (IsUndr .EQ. 1) Then
          ! BulkW = ...
          E_fac1      = max(c*cos(phi)+Prs_0(3)*sin(phi), max(c*cos(phi)-sigma_trac*sin(phi), 0d0))
          E_fac2      = (c*cos(phi)+p_ref*sin(phi))
          Eur         = Eur_ref*(E_fac1/E_fac2)**xm            
          G           = Eur/(2d0*(1d0+xNu_ur))      
          xNu_U   = 0.495d0
          Fac     = (1d0+xNu_U)/(1d0-2d0*xNu_U)-(1d0+xNu_ur)/(1d0-2d0*xNu_ur)
          Fac     = 2d0*G/3d0*Fac
          BulkW   = Fac
        End If
      End If  ! IDTask = 3, 6

    
      If (IDTask .EQ. 4) Then !Número de variables de estado
          
        nStat    = nStatV
        
      End If  ! IDTask = 4

    
      If (IDTask .Eq. 5) Then !Tipo de matriz
          
        NonSym   = 0  ! 1 para matrices D no simétricas
        iStrsDep = 1  ! 1 para matrices D que dependen de las tensiones
        iTang    = 0  ! 1 para matrices D tangentes
        iTimeDep = 0  ! 1 para matrices D que dependen del tiempo
        
      End If  ! IDTask = 5
    
      End ! MyModel1
____________________________________________________________________________________________________
 
 
HS_CAP.F90  (SUBROUTINE)
_____________________________________________________________________________________________________
!C***********************************************************************
      Subroutine HS_CAP_plast(iArea, Props, Prs_E, Prs, StVar0, StVar, ipl, IsUndr)
!C***********************************************************************
      Implicit Double Precision (A-H,O-Z)
!
! Rutina para el Returning Map en superficie HS
! Tensión de compresión positiva!
!
!        I/O  Tipo
! iArea   I    I  1 : Superficie f13_hs             Prs_E(1) > Prs_E(2) > Prs_E(3)
!                 2 : Eje entre f13_hs y f12_hs     Prs_E(1) > Prs_E(2) = Prs_E(3)
!                 3 : Eje entre f13_hs y f23_hs     Prs_E(1) = Prs_E(2) > Prs_E(3)
!      
! Local:                
!       f13_hs            : superficie HS para         Prs_E(1) > Prs_E(2) > Prs_E(3)
!       f12_hs            : superficie HS para         Prs_E(1) > Prs_E(2) = Prs_E(3)
!       f23_hs            : superficie HS para         Prs_E(1) = Prs_E(2) > Prs_E(3)
!       xLam13_hs         : multiplicador plástico f13_hs
!       xLam12_hs         : multiplicador plástico f12_hs
!       xLam23_hs         : multiplicador plástico f23_hs     
!       xLam_cap          : multiplicador plástico f_cap      
!       D                 : matriz elástica
!       df13_hsmdS(3)     : Primera derivada función de fluencia f13_hs/dS
!       df12_hsmdS(3)     : Primera derivada función de fluencia f12_hs/dS
!       df23_hsmdS(3)     : Primera derivada función de fluencia f23_hs/dS     
!       df_capmdS(3)      : Primera derivada función de fluencia f_cap/dS     
!       df13_hsmdq        : Primera derivada función de fluencia f13_hs/dq
!       df12_hsmdq        : Primera derivada función de fluencia f12_hs/dq
!       df23_hsmdq        : Primera derivada función de fluencia f23_hs/dq      
!       df_capmdq         : Primera derivada función de fluencia f_cap/dq      
!       dg13_hsmdS(3)     : Primera derivada potencial plástico g13_hs/dS
!       dg12_hsmdS(3)     : Primera derivada potencial plástico g12_hs/dS     
!       dg23_hsmdS(3)     : Primera derivada potencial plástico g23_hs/dS     
!       dg_capmdS(3)      : Primera derivada potencial plástico g_cap/dS           
!       d2g13_hsmdS2(3,3) : Segunda derivada potencial plástico g13_hs/dS2
!       d2g12_hsmdS2(3,3) : Segunda derivada potencial plástico g12_hs/dS2
!       d2g23_hsmdS2(3,3) : Segunda derivada potencial plástico g23_hs/dS2      
!       d2g_capmdS2(3,3)  : Segunda derivada potencial plástico g_cap/dS2            
!       d2g13_hsmdSdq(3)  : Derivada dS y dq del potencial plástico g13_hs       
!       d2g12_hsmdSdq(3)  : Derivada dS y dq del potencial plástico g12_hs       
!       d2g23_hsmdSdq(3)  : Derivada dS y dq del potencial plástico g23_hs       
!       d2g_capmdSdq(3)   : Derivada dS y dq del potencial plástico g_cap             
!       D_dgij            : Producto de D por dgij_hs/dS
!       D_dgij_2          : Producto de D por d2gij_hs/dS2
!       HH_dfij           : Producto de D por dfij_hs/dS
!       thetaij           : Operación dgij_hs/dS+d2gij_hs/dS2*xLamij_hs*xh_hs
!       D_thetaij         : Productor de D por thetaij
!       HH                : Inversa de la matriz Unity+xLamij_hs*D_gij_2   
      
      Dimension Prs_E(3), Prs(3), Prs_res1(3), Prs_res2(3), Prs_res3(3), Prs_0(3), xN10(3), xN20(3), xN30(3)
      Dimension Prs_2sup(3), Prs_3sup(3), r(3)
      Dimension Props(*)
      Dimension StVar0(*), StVar(*)      
      Dimension df13_hsmdS(3), df12_hsmdS(3), df23_hsmdS(3), df_capmdS(3) 
      Dimension dg13_hsmdS(3), dg12_hsmdS(3), dg23_hsmdS(3),  dg_capmdS(3)
      Dimension d2g13_hsmdS2(3,3), d2g12_hsmdS2(3,3), d2g23_hsmdS2(3,3), d2g_capmdS2(3,3)
      Dimension d2g13_hsmdSdq(3), d2g12_hsmdSdq(3), d2g23_hsmdSdq(3), d2g_capmdSdq(3)
      Dimension D_dg13(3), D_dg12(3), D_dg23(3), D_dgcap(3)
      Dimension D_dg13_2(3,3), D_dg12_2(3,3), D_dg23_2(3,3), D_dgcap_2(3,3) 
      Dimension HH_df13(3), HH_df12(3), HH_df23(3), HH_dfcap(3)
      Dimension theta13(3), theta12(3), theta23(3), thetacap(3)
      Dimension D_theta13(3), D_theta12(3), D_theta23(3), D_thetacap(3)
      Dimension Unity(3,3), D(3,3), HH(3,3), HH_inv(3,3), xMat(3,3), xMat_inv(3,3), f_b(3), d_lam(3)
      
      Call MZeroR(Prs_res1, 3)
      Call MZeroR(Prs_res2, 3)      
      Call MZeroR(Prs_res3, 3)      
      Call MZeroR(Prs_2sup, 3)
      Call MZeroR(Prs_3sup, 3)
            
      Call MZeroR(df13_hsmdS, 3) 
      Call MZeroR(dg13_hsmdS, 3)
      Call MZeroR(d2g13_hsmdS2, 9) 
      Call MZeroR(d2g13_hsmdSdq, 3)       

      Call MZeroR(df12_hsmdS, 3) 
      Call MZeroR(dg12_hsmdS, 3)
      Call MZeroR(d2g12_hsmdS2, 9) 
      Call MZeroR(d2g12_hsmdSdq, 3)       

      Call MZeroR(df23_hsmdS, 3) 
      Call MZeroR(dg23_hsmdS, 3)
      Call MZeroR(d2g23_hsmdS2, 9) 
      Call MZeroR(d2g23_hsmdSdq, 3)             
      
      Call MZeroR(df_capmdS, 3) 
      Call MZeroR(dg_capmdS, 3)
      Call MZeroR(d2g_capmdS2, 9) 
      Call MZeroR(d2g_capmdSdq, 3)             
      
      Call MZeroR(D_dg13_2, 9) 
      Call MZeroR(D_dg12_2, 9) 
      Call MZeroR(D_dg23_2, 9) 
      Call MZeroR(D_dgcap_2, 9) 
      
      Call MZeroR(D_dg13, 3) 
      Call MZeroR(D_dg12, 3) 
      Call MZeroR(D_dg23, 3) 
      Call MZeroR(D_dgcap, 3) 
             
      Call MZeroR(HH_df13, 3) 
      Call MZeroR(HH_df12, 3) 
      Call MZeroR(HH_df23, 3) 
      Call MZeroR(HH_dfcap, 3) 

      Call MZeroR(theta13, 3) 
      Call MZeroR(theta12, 3) 
      Call MZeroR(theta23, 3) 
      Call MZeroR(thetacap, 3) 

      Call MZeroR(D_theta13, 3) 
      Call MZeroR(D_theta12, 3) 
      Call MZeroR(D_theta23, 3) 
      Call MZeroR(D_thetacap, 3) 

      Call MZeroR(HH, 3) 
      Call MZeroR(HH_inv, 3) 
      Call MZeroR(xMat, 9)
      Call MZeroR(xMat_inv, 9)
      Call MZeroR(f_b, 3)      
      Call MZeroR(d_lam, 3)            
       
      ipl       = 5
      it        = 0
      itmax     = 20
      ttol      = 1d-9
      pi        = 3.1415926535d0
      Rad       = pi/180d0
      ideb      = -1 !para debugging      
      
      Call MZeroR(Unity, 9)
      Unity(1,1) = 1d0
      Unity(2,2) = 1d0
      Unity(3,3) = 1d0
      
      G         = Props(1) 
      xNu       = Props(2)     
      E50_ref   = Props(3) 
      Eoed_ref  = Props(4)
      Eur_ref   = Props(5)
      xm        = Props(6)
      c         = Props(7)
      phi       = Props(8)*Rad
      psi       = Props(9)*Rad
      xNu_ur    = Props(10)
      p_ref     = Props(11)
      xK0_nc    = Props(12)
      R_f       = Props(13)
      sigma_trac= Props(14)    
      
      xKs        = Eur_ref/((3d0*(1d0-2d0*xNu_ur))) 
      xKsmKc     = (Eur_ref/Eoed_ref)*(xK0_nc/((1d0+2d0*xK0_nc)*(1d0-2d0*xNu_ur)))
      xH         = xKs/(xKsmKc-1d0)
      
      gamma_p_E = StVar0(1)
      p_p_E     = (-1d0)*StVar0(2)
      
      dxLam13_hs = 0d0
      dxLam12_hs = 0d0
      dxLam23_hs = 0d0
      dxLam_cap  = 0d0

      xLam13_hs = 0d0
      xLam12_hs = 0d0
      xLam23_hs = 0d0
      xLam_cap  = 0d0
            
      !Construcción matriz D    
      iOpt        = 0
      Call PrnSig(iOpt, Sig0, xN10, xN20, xN30, S10, S20, S30, P0, Q0)
      Prs_0(1)    = S10*(-1d0)
      Prs_0(2)    = S20*(-1d0)
      Prs_0(3)    = S30*(-1d0)        
      Call Matriz3(Props, Prs_0, D, IsUndr)          
  
      Select Case (IArea)
      
      Case (1) ! 1 : Superficie f13_hs y f_cap (Sig1 > Sig2 > Sig3)         
         
         !STARTING POINT
         
         !Cálculo valor de f13_hs y f_cap en Prs_E         
         qf1         = max(c*cos(phi)/sin(phi)+Prs_E(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))!OJO, Prs(3) debe ser mayor que c*cos(phi)/sin(phi) para que E sean reales
         qf2         = (2d0*sin(phi)/(1d0-sin(phi)))
         qf          = qf1*qf2
         qa          = qf/R_f
         E_fac1      = max(c*cos(phi)+Prs_E(3)*sin(phi), c*cos(phi))
         E_fac2      = (c*cos(phi)+p_ref*sin(phi))
         E50         = E50_ref*(E_fac1/E_fac2)**xm
         Ei          = 2d0*E50/(2d0-R_f)  
         Eur         = Eur_ref*(E_fac1/E_fac2)**xm            
                 
         f13_hs_fac1 = ((Prs_E(1)-Prs_E(3))/(qa-(Prs_E(1)-Prs_E(3))))
         f13_hs_fac2 = 2d0*(Prs_E(1)-Prs_E(3))/Eur
         f13_hs      = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p_E
                   
         xlam_ast    = p_ref/Eoed_ref
         xkap_ast    = 3d0*p_ref*(1-2*xNu_ur)/Eur_ref
         alfa1       = ((1d0-xK0_nc)**2d0)/((1d0+2d0*xK0_nc)**2d0)
         alfa21      = (1d0-xK0_nc)*(1d0-2d0*xNu_ur)*(xlam_ast/xkap_ast-1d0)
         alfa221     = (1d0+2d0*xK0_nc)*(1d0-2d0*xNu_ur)*(xlam_ast/xkap_ast)
         alfa222     = (1d0-xK0_nc)*(1d0+xNu_ur)
         alfa        = 3d0*(alfa1+(alfa21/(alfa221-alfa222)))**0.5d0              
         delta       = (3d0+sin(phi))/(3d0-sin(phi))
         q_m         = Prs_E(1)+(delta-1d0)*Prs_E(2)-delta*Prs_E(3)
         p           = (Prs_E(1)+Prs_E(2)+Prs_E(3))/3d0          
 
         f_cap       = (q_m**2d0)/(alfa**2d0)+p**2d0-p_p_E**2d0  
         
         Call WriVal(ideb, 'qa', qa) 
         Call WriVal(ideb, 'qf', qf)    
         Call WriVal(ideb, 'Prs_E(1)-Prs_E(3)', Prs_E(1)-Prs_E(3))             
         Call WriVal(ideb, 'qa-(Prs_E(1)-Prs_E(3))', qa-(Prs_E(1)-Prs_E(3)))             
         
         !Cálculo tolerancia de f13_hs y f_cap
         f13_hs_TOL = ttol*abs(f13_hs)+ttol
         f_cap_TOL = ttol*abs(f_cap)+ttol

         Call WriVal(ideb, 'f13_hs', f13_hs)                                          
         Call WriVal(ideb, 'f_cap', f_cap)                                          
         Call WriVal(ideb, 'f13_hs_TOL', f13_hs_TOL)                                          
         Call WriVal(ideb, 'f_cap_TOL', f_cap_TOL)                                          
         
         !Calculo df13_hs/dS y df_cap/dS en Prs_E
         df13_hsmdS(1)    = ((2d0*(qa**2d0))/Ei)/((qa-(Prs_E(1)-Prs_E(3)))**2d0)-2d0/Eur
         df13_hsmdS(2)    = 0d0
         df13_hsmdS(3)    = -((2d0*(qa**2d0))/Ei)/((qa-(Prs_E(1)-Prs_E(3)))**2d0)+2d0/Eur             
         
         df_capmdS(1)    = 2d0*(q_m/(alfa**2d0)+(1d0/3d0)*p)
         df_capmdS(2)    = 2d0*(q_m*(delta-1d0)/(alfa**2d0)+(1d0/3d0)*p)
         df_capmdS(3)    = 2d0*(-q_m*delta/(alfa**2d0)+(1d0/3d0)*p)
         
         Call WriVec(ideb, 'df13_hsmdS', df13_hsmdS, 3)                                          
         Call WriVec(ideb, 'df_capmdS', df_capmdS, 3)                                                   
         
         !Cálculo df13_hs/dq en Prs_E (q1=gamma_p, q2=p_p)
         df13_hsmdq      = -1d0
         df_capmdq       = -2d0*p_p_E   
        
         !Cálculo h en Prs_E (dq1=dgamma_p=dlam1*h1, dq2=dp_p=dlam2*h2)
         xh_hs           = 1d0
         xh_cap          = 2d0*xH*p*(max(Prs_E(3)+c*cos(phi)/sin(phi),c*cos(phi)/sin(phi))/(p_ref+c*cos(phi)/sin(phi)))**xm
         
         !Cálculo dg13_hs/dS y dg_cap/dS en Prs_E     
         sin_phi_m        = (Prs_E(1)-Prs_E(3))/(Prs_E(1)+Prs_E(3)-2d0*c*cos(phi)/sin(phi))
         sin_phi_cv       = (sin(phi)-sin(psi))/(1d0-(sin(phi))*(sin(psi)))
           If (sin_phi_m .LT. (3d0/4d0)*sin(phi)) Then
                sin_psi_m = 0d0
           ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .GT. 0d0) Then
                sin_psi_m = max((sin_phi_m-sin_phi_cv)/(1d0-sin_phi_m*sin_phi_cv),0d0)
           ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .LE. 0d0) Then
                psi_m = psi
                sin_psi_m = sin(psi_m)
           ElseIf (phi .EQ. 0d0) Then
                sin_psi_m = 0d0
           End If
         
         dg13_hsmdS(1)   = (1d0/2d0)*(1d0-sin_psi_m)  
         dg13_hsmdS(2)   = 0d0
         dg13_hsmdS(3)   = -(1d0/2d0)*(1d0+sin_psi_m)         
 
         dg_capmdS(1)   = df_capmdS(1) 
         dg_capmdS(2)   = df_capmdS(2) 
         dg_capmdS(3)   = df_capmdS(3) 

         Call WriVec(ideb, 'dg13_hsmdS', dg13_hsmdS, 3)                                          
         Call WriVec(ideb, 'dg_capmdS', dg_capmdS, 3)                                         
         
         !Cálculo d2g13_hs/dS2 y d2gcap/dS2 en Prs_E   
         Call MZeroR(d2g13_hsmdS2, 9) 

         d2g_capmdS2(1,1) = 2d0*(1d0/(alfa**2d0)+1d0/9d0)
         d2g_capmdS2(1,2) = 2d0*((delta-1d0)/(alfa**2d0)+1d0/9d0)
         d2g_capmdS2(1,3) = -2d0*(delta/(alfa**2d0)-1d0/9d0)
         d2g_capmdS2(2,1) = d2g_capmdS2(1,2)    
         d2g_capmdS2(2,2) = 2d0*(((delta-1d0)**2d0)/(alfa**2d0)+1d0/9d0)   
         d2g_capmdS2(2,3) = -2d0*(((delta-1d0)*delta)/(alfa**2d0)-1d0/9d0)   
         d2g_capmdS2(3,1) = d2g_capmdS2(1,3)   
         d2g_capmdS2(3,2) = d2g_capmdS2(2,3)   
         d2g_capmdS2(3,3) = 2d0*((delta**2d0)/(alfa**2d0)+1d0/9d0)   
         
         !Cálculo d2g13_hs/dSdq en Prs_E
         Call MZeroR(d2g13_hsmdSdq, 3) 
         Call MZeroR(d2g_capmdSdq, 3)          
         
        !Cálculo de los multiplicadores plásticos 
         Call MatMat(D, 3, dg13_hsmdS, 3, 3, 1, 3, D_dg13, 3) 
         Call DInProd(df13_hsmdS, D_dg13, 3, df13_D_dg13)
         xOmega11 = df13_D_dg13-df13_hsmdq*xh_hs        
         Call MatMat(D, 3, dg_capmdS, 3, 3, 1, 3, D_dgcap, 3) 
         Call DInProd(df13_hsmdS, D_dgcap, 3, df13_D_dgcap)
         xOmega12 = df13_D_dgcap      
         Call DInProd(df_capmdS, D_dg13,3, dfcap_D_dg13)
         xOmega21 = dfcap_D_dg13                     
         Call DInProd(df_capmdS, D_dgcap,3, dfcap_D_dgcap)
         xOmega22 = dfcap_D_dgcap-df_capmdq*xh_cap        
         
!         Call WriVal(ideb, 'df13_D_dg13', df13_D_dg13)                                 
!         Call WriVal(ideb, 'df13_hsmdq', df13_hsmdq)             
!         Call WriVal(ideb, 'xh_hs', xh_hs)             
!         Call WriVal(ideb, 'xOmega11 = df13_D_dg13-df13_hsmdq*xh_hs', df13_D_dg13-df13_hsmdq*xh_hs)                                          
!         Call WriVal(ideb, 'xOmega12 = df13_D_dgcap', df13_D_dgcap)                                              
!         Call WriVal(ideb, 'xOmega21 = dfcap_D_dg13p', dfcap_D_dg13)  
!         Call WriVal(ideb, 'dfcap_D_dgcap', dfcap_D_dgcap)                                 
!         Call WriVal(ideb, 'df_capmdq', df_capmdq)             
!         Call WriVal(ideb, 'xh_cap', xh_cap)             
!         Call WriVal(ideb, 'xOmega22 = dfcap_D_dgcap-df_capmdq*xh_cap ', dfcap_D_dgcap-df_capmdq*xh_cap)                                     
         
         !Resuolución sistema y cálculo dxLam13_hs y dxLam_cap
         dxLam13_hs = (f13_hs*xOmega22-f_cap*xOmega12)/(xOmega11*xOmega22-xOmega12*xOmega21)
         dxLam_cap = (f_cap*xOmega11-f13_hs*xOmega21)/(xOmega11*xOmega22-xOmega12*xOmega21)             
         
         xLam13_hs = dxLam13_hs         
         xLam_cap = dxLam_cap         

         Call WriVal(ideb, 'dxLam13_hs', dxLam13_hs)                                 
         Call WriVal(ideb, 'dxLam_cap', dxLam_cap)                                              
         Call WriVal(ideb, 'xLam13_hs', xLam13_hs)                                 
         Call WriVal(ideb, 'xLam_cap', xLam_cap)                                              
                  
         !Actualizo variables
         Call AddVec(D_dg13, D_dgcap, -xLam13_hs, -xLam_cap, 3, Prs_2sup)         
         Call AddVec(Prs_E, Prs_2sup, 1d0, 1d0, 3, Prs)
         gamma_p = gamma_p_E + xLam13_hs*xh_hs
         p_p = p_p_E + xLam_cap*xh_cap

         !Cálculo valor de f13_hs y f_cap con la nueva Prs                  
         qf1         = max(c*cos(phi)/sin(phi)+Prs(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))  !OJO, Prs(3) debe ser mayor que c*cos(phi)/sin(phi) para que E sean reales
         qf2         = (2d0*sin(phi)/(1d0-sin(phi)))
         qf          = qf1*qf2
         qa          = qf/R_f
         E_fac1      = max(c*cos(phi)+Prs(3)*sin(phi), c*cos(phi))
         E_fac2      = (c*cos(phi)+p_ref*sin(phi))
         E50         = E50_ref*(E_fac1/E_fac2)**xm
         Ei          = 2d0*E50/(2d0-R_f)
         Eur         = Eur_ref*(E_fac1/E_fac2)**xm          

         f13_hs_fac1 = ((Prs(1)-Prs(3))/(qa-(Prs(1)-Prs(3))))
         f13_hs_fac2 = 2d0*(Prs(1)-Prs(3))/Eur
         f13_hs      = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p   
                                  
         q_m         = Prs(1)+(delta-1d0)*Prs(2)-delta*Prs(3)
         p           = (Prs(1)+Prs(2)+Prs(3))/3d0          
         
         f_cap       = (q_m**2d0)/(alfa**2d0)+p**2d0-p_p**2d0            
         
         !EULER BACKWARD
         
         Do while ((abs(f13_hs) .GT. f13_hs_TOL .OR. &
                   abs(f_cap) .GT. f_cap_TOL) .AND. &
                   it .LE. itmax)

             !Construcción matriz D
             !Call Matriz3(Props, Prs, D, IsUndr)                          
             
             !Calculo df13_hs/dS y df_cap/dS en Prs
         
             Call WriVal(ideb, 'qa', qa) 
             Call WriVal(ideb, 'qf', qf)    
             Call WriVal(ideb, 'Prs(1)-Prs(3)', Prs(1)-Prs(3))             
             Call WriVal(ideb, 'qa-(Prs(1)-Prs(3))', qa-(Prs(1)-Prs(3)))             
            
             df13_hsmdS(1)    = ((2d0*(qa**2d0))/Ei)*(1/((qa-(Prs(1)-Prs(3)))**2d0))-2d0/Eur
             df13_hsmdS(2)    = 0d0
             df13_hsmdS(3)    = -((2d0*(qa**2d0))/Ei)*(1/((qa-(Prs(1)-Prs(3)))**2d0))+2d0/Eur    
                 
             df_capmdS(1)    = 2d0*(q_m/(alfa**2d0)+(1d0/3d0)*p)
             df_capmdS(2)    = 2d0*(q_m*(delta-1d0)/(alfa**2d0)+(1d0/3d0)*p)
             df_capmdS(3)    = 2d0*(-q_m*delta/(alfa**2d0)+(1d0/3d0)*p)
          
             Call WriVec(ideb, 'df13_hsmdS', df13_hsmdS, 3)             
             Call WriVec(ideb, 'df_capmdS', df_capmdS, 3)                    
                         
             !Cálculo df13_hs/dq en Prs (q1=gamma_p, q2=p_p)
             df13_hsmdq      = -1d0
             df_capmdq       = -2d0*p_p   
             
             Call WriVal(ideb, 'df_capmdq', df_capmdq)                    
             
             !Cálculo h en Prs (dq=dlam1*h1 y dq2=dlam2*h2)
             xh_hs           = 1d0
             xh_cap          = 2d0*xH*p*(max(Prs(3)+c*cos(phi)/sin(phi), c*cos(phi)/sin(phi))/(p_ref+c*cos(phi)/sin(phi)))**xm          

             Call WriVal(ideb, 'xh_cap', xh_cap)                                 
             
             !Cálculo dg13_hs/dS y dg_cap/dS en Prs     
             sin_phi_m        = (Prs(1)-Prs(3))/(Prs(1)+Prs(3)-2d0*c*cos(phi)/sin(phi))
             sin_phi_cv       = (sin(phi)-sin(psi))/(1d0-(sin(phi))*(sin(psi)))
             If (sin_phi_m .LT. (3d0/4d0)*sin(phi)) Then
                 sin_psi_m = 0d0
             ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .GT. 0d0) Then
                 sin_psi_m = max((sin_phi_m-sin_phi_cv)/(1d0-sin_phi_m*sin_phi_cv),0d0)
             ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .LE. 0d0) Then
                 psi_m = psi
                 sin_psi_m = sin(psi_m)
             ElseIf (phi .EQ. 0d0) Then
                 sin_psi_m = 0d0
             End If
        
             dg13_hsmdS(1)   = (1d0/2d0)*(1d0-sin_psi_m)  
             dg13_hsmdS(2)   = 0d0
             dg13_hsmdS(3)   = -(1d0/2d0)*(1d0+sin_psi_m)         
              
             dg_capmdS(1)   = df_capmdS(1) 
             dg_capmdS(2)   = df_capmdS(2) 
             dg_capmdS(3)   = df_capmdS(3) 
                         
             Call WriVec(ideb, 'dg13_hsmdS', dg13_hsmdS, 3)                                          
             Call WriVec(ideb, 'dg_capmdS', dg_capmdS, 3)                                                   
             
             !Cálculo d2g13_hs/dS2 y d2g_cap/dS2 en Prs   
             Call MZeroR(d2g13_hsmdS2, 9) 
             
             d2g_capmdS2(1,1) = 2d0*(1d0/(alfa**2d0)+1d0/9d0)
             d2g_capmdS2(1,2) = 2d0*((delta-1d0)/(alfa**2d0)+1d0/9d0)
             d2g_capmdS2(1,3) = -2d0*(delta/(alfa**2d0)-1d0/9d0)
             d2g_capmdS2(2,1) = d2g_capmdS2(1,2)    
             d2g_capmdS2(2,2) = 2d0*(((delta-1d0)**2d0)/(alfa**2d0)+1d0/9d0)   
             d2g_capmdS2(2,3) = -2d0*(((delta-1d0)*delta)/(alfa**2d0)-1d0/9d0)   
             d2g_capmdS2(3,1) = d2g_capmdS2(1,3)   
             d2g_capmdS2(3,2) = d2g_capmdS2(2,3)   
             d2g_capmdS2(3,3) = 2d0*((delta**2d0)/(alfa**2d0)+1d0/9d0)   
          
             !Cálculo d2g13_hs/dSdq en Prs
             Call MZeroR(d2g13_hsmdSdq, 3) 
             Call MZeroR(d2g_capmdSdq, 3) 
         
             !Cálculo del multiplicador plástico 
             Call MatMat(D, 3, d2g13_hsmdS2, 3, 3, 3, 3, D_dg13_2, 3)
             Call MatMat(D, 3, d2g_capmdS2, 3, 3, 3, 3, D_dgcap_2, 3)             
             HH_inv = Unity + xLam13_hs*D_dg13_2 + xLam_cap*D_dgcap_2
             Call MatInvPiv(HH_inv, HH, 3)
             Call AddVec(dg13_hsmdS, d2g13_hsmdSdq, 1d0, xLam13_hs*xh_hs, 3, theta13) 
             Call AddVec(dg_capmdS, d2g_capmdSdq, 1d0, xLam_cap*xh_cap, 3, thetacap) 
             Call MatMat(HH, 3, df13_hsmdS, 3, 3, 1, 3, HH_df13,3)
             Call MatMat(HH, 3, df_capmdS, 3, 3, 1, 3, HH_dfcap,3)             
             Call MatMat(D, 3, theta13, 3, 3, 1, 3, D_theta13,3)
             Call MatMat(D, 3, thetacap, 3, 3, 1, 3, D_thetacap,3)             
             Call DInProd(D_theta13, HH_df13, 3, D_theta13_HH_df13) 
             Call DInProd(D_thetacap, HH_df13, 3, D_thetacap_HH_df13)          
             Call DInProd(D_theta13, HH_dfcap, 3, D_theta13_HH_dfcap)                           
             Call DInProd(D_thetacap, HH_dfcap, 3, D_thetacap_HH_dfcap)    
             xOmega11  = D_theta13_HH_df13-df13_hsmdq*xh_hs                     
             xOmega12  = D_thetacap_HH_df13                                 
             xOmega21  = D_theta13_HH_dfcap                              
             xOmega22  = D_thetacap_HH_dfcap-df_capmdq*xh_cap    

             !Call WriMat(1, 'HH', HH, 3, 3, 3)  
             Call MatMat(D, 3, dg13_hsmdS, 3, 3, 1, 3, D_dg13, 3)       
             Call MatMat(D, 3, dg_capmdS, 3, 3, 1, 3, D_dgcap, 3)                      
             Call AddVec(D_dg13, D_dgcap, -xLam13_hs, -xLam_cap, 3, Prs_res1)
             Call AddVec(Prs_E, Prs_res1, 1d0, 1d0, 3, Prs_res2)
             Call AddVec(Prs, Prs_res2, 1d0, -1d0, 3, r)         
             Call DInProd(HH_df13, r, 3, r_HH_df13)                           
             Call DInProd(HH_dfcap, r, 3, r_HH_dfcap)             
             xOmegaf13 = f13_hs - r_HH_df13
             xOmegafcap = f_cap - r_HH_dfcap             
             
             !Resolución sistema y cálculo dxLam13_hs y dxLam_cap            
             dxLam13_hs = (xOmegaf13*xOmega22-xOmegafcap*xOmega12)/(xOmega11*xOmega22-xOmega12*xOmega21)
             dxLam_cap  = (xOmegafcap*xOmega11-xOmegaf13*xOmega21)/(xOmega11*xOmega22-xOmega12*xOmega21)              
     
             xLam13_hs = xLam13_hs + dxLam13_hs
             xLam_cap  = xLam_cap + dxLam_cap
             
             Call WriVal(ideb, 'dxLam13_hs', dxLam13_hs)                                 
             Call WriVal(ideb, 'dxLam_cap', dxLam_cap)                                              
             Call WriVal(ideb, 'xLam13_hs', xLam13_hs)                                 
             Call WriVal(ideb, 'xLam_cap', xLam_cap)                                              
             
             !Actualización variables
             Call AddVec(D_dg13, D_dgcap, -xLam13_hs, -xLam_cap, 3, Prs_2sup)         
             Call AddVec(Prs_E, Prs_2sup, 1d0, 1d0, 3, Prs)
             gamma_p = gamma_p_E + xLam13_hs*xh_hs
             p_p = p_p_E + xLam_cap*xh_cap

             Call WriVal(ideb, 'gamma_p', gamma_p)                                 
             Call WriVal(ideb, 'p_p', p_p)                                              
             
             Call WriVal(ideb, 'qa', qa)    
             Call WriVal(ideb, 'Prs(1)-Prs(3)', Prs(1)-Prs(3))             
             Call WriVal(ideb, 'qa-(Prs(1)-Prs(3))', qa-(Prs(1)-Prs(3)))                   
             Call WriVal(ideb, 'Ei', Ei)     
             Call WriVal(ideb, 'Eur', Eur)                  
             Call WriVal(ideb, '(qa-(Prs(1)-Prs(3)))**2d0', (qa-(Prs(1)-Prs(3)))**2d0)     
             
             Call WriVec(ideb, 'Prs durante iteración', Prs, 3)                  
             
             !Cálculo valor de f13_hs en (k+1) 
             qf1         = max(c*cos(phi)/sin(phi)+Prs(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))!OJO, Prs(3) debe ser mayor que c*cos(phi)/sin(phi) para que E sean reales
             qf2         = (2d0*sin(phi)/(1d0-sin(phi)))
             qf          = qf1*qf2
             qa          = qf/R_f
             E_fac1      = max(c*cos(phi)+Prs(3)*sin(phi), c*cos(phi))
             E_fac2      = (c*cos(phi)+p_ref*sin(phi))
             E50         = E50_ref*(E_fac1/E_fac2)**xm
             Ei          = 2d0*E50/(2d0-R_f)
             Eur         = Eur_ref*(E_fac1/E_fac2)**xm          
             
             f13_hs_fac1 = ((Prs(1)-Prs(3))/(qa-(Prs(1)-Prs(3))))
             f13_hs_fac2 = 2d0*(Prs(1)-Prs(3))/Eur
             f13_hs      = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p   
                                       
             q_m         = Prs(1)+(delta-1d0)*Prs(2)-delta*Prs(3)
             p           = (Prs(1)+Prs(2)+Prs(3))/3d0          
         
             f_cap       = (q_m**2d0)/(alfa**2d0)+p**2d0-p_p**2d0                   

            Call WriIvl(ideb, 'it', it)                
             
             !Actualizo iteración
             it = it+1
  
         End do
            
         !Actualizo variables de estado
         StVar(1) = gamma_p
         StVar(2) = min((-1d0)*p_p, (-1d0)*p_p_E)

        Call WriVal(ideb, 'xLam13_hs', xLam13_hs)                
        Call WriVal(ideb, 'xLam_cap', xLam_cap)                
        
         If (xLam13_hs .GT. (-1d0)*ttol .AND. &
             xLam_cap .GT. (-1d0)*ttol .AND. &
             it .LE. itmax) Then
             write(1,*)'Return mapping Ok!'
         Else
             write(1,*)'Numerical problem'
         End If
         
      
      Case (2) !  2 : Eje entre f13_hs, f12_hs y f_cap (Prs(1) > Prs(2) = Prs(3))

         !STARTING POINT
         
         !Cálculo el valor de f13_hs, f12_hs y f_cap en Prs_E
         qf1         = max(c*cos(phi)/sin(phi)+Prs_E(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))!OJO, Prs(3) debe ser mayor que c*cos(phi)/sin(phi) para que E sean reales
         qf2         = (2d0*sin(phi)/(1d0-sin(phi)))
         qf          = qf1*qf2
         qa          = qf/R_f
         E_fac1      = max(c*cos(phi)+Prs_E(3)*sin(phi), c*cos(phi))
         E_fac2      = (c*cos(phi)+p_ref*sin(phi))
         E50         = E50_ref*(E_fac1/E_fac2)**xm
         Ei          = 2d0*E50/(2d0-R_f)
         Eur         = Eur_ref*(E_fac1/E_fac2)**xm          
         
         f13_hs_fac1 = ((Prs_E(1)-Prs_E(3))/(qa-(Prs_E(1)-Prs_E(3))))
         f13_hs_fac2 = 2d0*(Prs_E(1)-Prs_E(3))/Eur
         f13_hs      = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p_E

         f12_hs_fac1 = ((Prs_E(1)-Prs_E(2))/(qa-(Prs_E(1)-Prs_E(2))))
         f12_hs_fac2 = 2d0*(Prs_E(1)-Prs_E(2))/Eur
         f12_hs      = (2d0*qa/Ei)*f12_hs_fac1-f12_hs_fac2-gamma_p_E
         
         xlam_ast    = p_ref/Eoed_ref
         xkap_ast    = 3d0*p_ref*(1-2*xNu_ur)/Eur_ref
         alfa1       = ((1d0-xK0_nc)**2d0)/((1d0+2d0*xK0_nc)**2d0)
         alfa21      = (1d0-xK0_nc)*(1d0-2d0*xNu_ur)*(xlam_ast/xkap_ast-1d0)
         alfa221     = (1d0+2d0*xK0_nc)*(1d0-2d0*xNu_ur)*(xlam_ast/xkap_ast)
         alfa222     = (1d0-xK0_nc)*(1d0+xNu_ur)
         alfa        = 3d0*(alfa1+(alfa21/(alfa221-alfa222)))**0.5d0             
         delta       = (3d0+sin(phi))/(3d0-sin(phi))
         q_m         = Prs_E(1)+(delta-1d0)*Prs_E(2)-delta*Prs_E(3)
         p           = (Prs_E(1)+Prs_E(2)+Prs_E(3))/3d0          

         f_cap       = (q_m**2d0)/(alfa**2d0)+p**2d0-p_p_E**2d0  

         Call WriVal(ideb, 'qa', qa) 
         Call WriVal(ideb, 'qf', qf)    
         Call WriVal(ideb, 'Prs_E(1)-Prs_E(3)', Prs_E(1)-Prs_E(3))             
         Call WriVal(ideb, 'qa-(Prs_E(1)-Prs_E(3))', qa-(Prs_E(1)-Prs_E(3)))             
         
         !Cálculo tolerancia f13_hs, f12_hs y f_cap en Prs_E
         f13_hs_TOL = ttol*abs(f13_hs)+ttol
         f12_hs_TOL = ttol*abs(f12_hs)+ttol
         f_cap_TOL = ttol*abs(f_cap)+ttol
         
         !Calculo df13_hs/dS, df12_hs/dS y df_cap/dS en Prs_E

         df13_hsmdS(1)    = ((2d0*(qa**2d0))/Ei)/((qa-(Prs_E(1)-Prs_E(3)))**2d0)-2d0/Eur
         df13_hsmdS(2)    = 0d0
         df13_hsmdS(3)    = -((2d0*(qa**2d0))/Ei)/((qa-(Prs_E(1)-Prs_E(3)))**2d0)+2d0/Eur 
         
         df12_hsmdS(1)    = ((2d0*(qa**2d0))/Ei)/((qa-(Prs_E(1)-Prs_E(2)))**2d0)-2d0/Eur
         df12_hsmdS(2)    = -((2d0*(qa**2d0))/Ei)/((qa-(Prs_E(1)-Prs_E(2)))**2d0)+2d0/Eur
         df12_hsmdS(3)    = 0d0  
                  
         df_capmdS(1)    = 2d0*(q_m/(alfa**2d0)+(1d0/3d0)*p)
         df_capmdS(2)    = 2d0*(q_m*(delta-1d0)/(alfa**2d0)+(1d0/3d0)*p)
         df_capmdS(3)    = 2d0*(-q_m*delta/(alfa**2d0)+(1d0/3d0)*p)         

         Call WriVec(ideb, 'df13_hsmdS', df13_hsmdS, 3)                                          
         Call WriVec(ideb, 'df12_hsmdS', df12_hsmdS, 3)                
         Call WriVec(ideb, 'df_capmdS', df_capmdS, 3)         
                  
         !Cálculo df13_hs/dq, df12_hs/dq y df_cap/dS en Prs_E (q1=gamma_p  -para f13_hs y f12_hs-, q2=p_p)
         df13_hsmdq      = -1d0
         df12_hsmdq      = -1d0
         df_capmdq       = -2d0*p_p_E            
         
         Call WriVal(ideb, 'df13_hsmdq', df13_hsmdq)                    
         Call WriVal(ideb, 'df12_hsmdq', df12_hsmdq)                    
         Call WriVal(ideb, 'df_capmdq', df_capmdq)                    
         
         !Cálculo h en Prs_E (dq=dlam1*h1 -para f13_hs y f12_hs- y dq2=dlam2*h2)
         xh_hs           = 1d0
         xh_cap          = 2d0*xH*p*(max(Prs_E(3)+c*cos(phi)/sin(phi), c*cos(phi)/sin(phi))/(p_ref+c*cos(phi)/sin(phi)))**xm     

         Call WriVal(ideb, 'xh_hs', xh_hs)                    
         Call WriVal(ideb, 'xh_cap', xh_cap)                    
         
         !Cálculo dg13_hs/dS, dg12_hs/dS y dg_cap/dS en Prs_E
         sin_phi_m        = (Prs_E(1)-Prs_E(3))/(Prs_E(1)+Prs_E(3)-2d0*c*cos(phi)/sin(phi))
         sin_phi_cv       = (sin(phi)-sin(psi))/(1d0-(sin(phi))*(sin(psi)))
           If (sin_phi_m .LT. (3d0/4d0)*sin(phi)) Then
                sin_psi_m = 0d0
           ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .GT. 0d0) Then
                sin_psi_m = max((sin_phi_m-sin_phi_cv)/(1d0-sin_phi_m*sin_phi_cv),0d0)
           ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .LE. 0d0) Then
                psi_m = psi
                sin_psi_m = sin(psi_m)
           ElseIf (phi .EQ. 0d0) Then
                sin_psi_m = 0d0
           End If
           
         dg13_hsmdS(1)   = (1d0/2d0)*(1d0-sin_psi_m)  
         dg13_hsmdS(2)   = 0d0
         dg13_hsmdS(3)   = -(1d0/2d0)*(1d0+sin_psi_m)         

         dg12_hsmdS(1)   = (1d0/2d0)*(1d0-sin_psi_m)  
         dg12_hsmdS(2)   = -(1d0/2d0)*(1d0+sin_psi_m)
         dg12_hsmdS(3)   = 0d0         

         dg_capmdS(1)   = df_capmdS(1) 
         dg_capmdS(2)   = df_capmdS(2) 
         dg_capmdS(3)   = df_capmdS(3) 

         Call WriVec(ideb, 'dg13_hsmdS', dg13_hsmdS, 3)                                          
         Call WriVec(ideb, 'dg12_hsmdS', dg12_hsmdS, 3)                                      
         Call WriVec(ideb, 'dg_capmdS', dg_capmdS, 3)                                         
         
         !Cálculo d2g13_hs/dS2, d2g12_hs/dS2 y d2g_cap/dS2 en Prs_E
         Call MZeroR(d2g13_hsmdS2, 9) 
         Call MZeroR(d2g12_hsmdS2, 9)          
         
         d2g_capmdS2(1,1) = 2d0*(1d0/(alfa**2d0)+1d0/9d0)
         d2g_capmdS2(1,2) = 2d0*((delta-1d0)/(alfa**2d0)+1d0/9d0)
         d2g_capmdS2(1,3) = -2d0*(delta/(alfa**2d0)-1d0/9d0)
         d2g_capmdS2(2,1) = d2g_capmdS2(1,2)    
         d2g_capmdS2(2,2) = 2d0*(((delta-1d0)**2d0)/(alfa**2d0)+1d0/9d0)   
         d2g_capmdS2(2,3) = -2d0*(((delta-1d0)*delta)/(alfa**2d0)-1d0/9d0)   
         d2g_capmdS2(3,1) = d2g_capmdS2(1,3)   
         d2g_capmdS2(3,2) = d2g_capmdS2(2,3)   
         d2g_capmdS2(3,3) = 2d0*((delta**2d0)/(alfa**2d0)+1d0/9d0)   
         
         !Cálculo d2g13_hs/dSdq, d2g12_hs/dSdq y d2g_cap/dSdq en Prs_E
         Call MZeroR(d2g13_hsmdSdq, 3) 
         Call MZeroR(d2g12_hsmdSdq, 3)  
         Call MZeroR(d2g_capmdSdq, 3)          
         
         !Cálculo de los multiplicadores plásticos 
         Call MatMat(D, 3, dg13_hsmdS, 3, 3, 1, 3, D_dg13, 3) 
         Call DInProd(df13_hsmdS, D_dg13,3, df13_D_dg13)
         xOmega11 = df13_D_dg13-df13_hsmdq*xh_hs        
         Call MatMat(D,3,dg12_hsmdS,3,3,1,3,D_dg12,3) 
         Call DInProd(df13_hsmdS, D_dg12,3, df13_D_dg12)
         xOmega12 = df13_D_dg12     
         Call MatMat(D,3,dg_capmdS,3,3,1,3,D_dgcap,3) 
         Call DInProd(df13_hsmdS, D_dgcap,3, df13_D_dgcap)
         xOmega13 = df13_D_dgcap              
         Call DInProd(df12_hsmdS, D_dg13,3, df12_D_dg13)
         xOmega21 = df12_D_dg13               
         Call DInProd(df12_hsmdS, D_dg12,3, df12_D_dg12)
         xOmega22 = df12_D_dg12-df12_hsmdq*xh_hs 
         Call DInProd(df12_hsmdS, D_dgcap,3, df12_D_dgcap)
         xOmega23 = df12_D_dgcap
         Call DInProd(df_capmdS, D_dg13,3, dfcap_D_dg13)
         xOmega31 = dfcap_D_dg13
         Call DInProd(df_capmdS, D_dg12,3, dfcap_D_dg12)
         xOmega32 = dfcap_D_dg12
         Call DInProd(df_capmdS, D_dgcap,3, dfcap_D_dgcap)
         xOmega33 = dfcap_D_dgcap-df_capmdq*xh_cap          
                         
         !Resolución sistema y calculo dxLam13_hs, dxLam12_hs y dxLam_cap
         xMat(1,1) = xOmega11
         xMat(1,2) = xOmega12
         xMat(1,3) = xOmega13         
         xMat(2,1) = xOmega21         
         xMat(2,2) = xOmega22         
         xMat(2,3) = xOmega23         
         xMat(3,1) = xOmega31         
         xMat(3,2) = xOmega32         
         xMat(3,3) = xOmega33         
         f_b(1) = f13_hs
         f_b(2) = f12_hs         
         f_b(3) = f_cap         
         
         Call MatInvPiv(xMat, xMat_inv, 3)
         Call MatMat(xMat_inv, 3, f_b, 3, 3, 1, 3, d_lam, 3) 
    
         dxLam13_hs = d_lam(1)
         dxLam12_hs = d_lam(2)
         dxLam_cap  = d_lam(3)
         
         xLam13_hs = dxLam13_hs         
         xLam12_hs = dxLam12_hs         
         xLam_cap  = dxLam_cap                  

         Call WriVal(ideb, 'dxLam13_hs', dxLam13_hs)                                 
         Call WriVal(ideb, 'dxLam12_hs', dxLam12_hs)                                          
         Call WriVal(ideb, 'dxLam_cap', dxLam_cap)                                              
         Call WriVal(ideb, 'xLam13_hs', xLam13_hs)                                 
         Call WriVal(ideb, 'xLam12_hs', xLam12_hs)                                          
         Call WriVal(ideb, 'xLam_cap', xLam_cap)                                              
                           
         !Actualización variables
         Call AddVec(D_dg13, D_dg12, -xLam13_hs, -xLam12_hs, 3, Prs_2sup) 
         Call AddVec(Prs_2sup, D_dgcap, 1d0, -xLam_cap, 3, Prs_3sup)          
         Call AddVec(Prs_E, Prs_3sup, 1d0, 1d0, 3, Prs)
         gamma_p = gamma_p_E + xLam13_hs*xh_hs
         p_p = p_p_E + xLam_cap*xh_cap

         Call WriVal(ideb, 'gamma_p', gamma_p)                                 
         Call WriVal(ideb, 'p_p', p_p)                                              
         
         !Cálculo valor de f13_hs, f12_hs y f_cap con la nueva Prs         
         qf1         = max(c*cos(phi)/sin(phi)+Prs(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))!OJO, Prs(3) debe ser mayor que c*cos(phi)/sin(phi) para que E sean reales
         qf2         = (2d0*sin(phi)/(1d0-sin(phi)))
         qf          = qf1*qf2
         qa          = qf/R_f
         E_fac1      = (c*cos(phi)+Prs(3)*sin(phi))
         E_fac2      = max(c*cos(phi)+p_ref*sin(phi), c*cos(phi))
         E50         = E50_ref*(E_fac1/E_fac2)**xm
         Ei          = 2d0*E50/(2d0-R_f)
         Eur         = Eur_ref*(E_fac1/E_fac2)**xm          

         f13_hs_fac1 = ((Prs(1)-Prs(3))/(qa-(Prs(1)-Prs(3))))
         f13_hs_fac2 = 2d0*(Prs(1)-Prs(3))/Eur
         f13_hs      = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p

         f12_hs_fac1 = ((Prs(1)-Prs(2))/(qa-(Prs(1)-Prs(2))))
         f12_hs_fac2 = 2d0*(Prs(1)-Prs(2))/Eur
         f12_hs      = (2d0*qa/Ei)*f12_hs_fac1-f12_hs_fac2-gamma_p
         
         q_m         = Prs(1)+(delta-1d0)*Prs(2)-delta*Prs(3)
         p           = (Prs(1)+Prs(2)+Prs(3))/3d0          
         
         f_cap       = (q_m**2d0)/(alfa**2d0)+p**2d0-p_p**2d0     
         
         !EULER BACKWARD
         
         Do while ((abs(f13_hs) .GT. f13_hs_TOL .OR. &
                   abs(f12_hs) .GT. f12_hs_TOL .OR. &
                   abs(f_cap) .GT. f_cap_TOL) .AND. &
                   it .LE. itmax)

             !Construcción matriz D
             !Call Matriz3(Props, Prs, D, IsUndr)                          

             Call WriVal(ideb, 'qa', qa) 
             Call WriVal(ideb, 'qf', qf)    
             Call WriVal(ideb, 'Prs(1)-Prs(3)', Prs(1)-Prs(3))             
             Call WriVal(ideb, 'qa-(Prs(1)-Prs(3))', qa-(Prs(1)-Prs(3)))             
            
             !Cálculo df13_hs/dS, df12_hs/dS y df_cap/dS en Prs
             df13_hsmdS(1)    = ((2d0*(qa**2d0))/Ei)/((qa-(Prs(1)-Prs(3)))**2d0)-2d0/Eur
             df13_hsmdS(2)    = 0d0
             df13_hsmdS(3)    = -((2d0*(qa**2d0))/Ei)/((qa-(Prs(1)-Prs(3)))**2d0)+2d0/Eur 
             
             df12_hsmdS(1)    = ((2d0*(qa**2d0))/Ei)/((qa-(Prs(1)-Prs(2)))**2d0)-2d0/Eur
             df12_hsmdS(2)    = -((2d0*(qa**2d0))/Ei)/((qa-(Prs(1)-Prs(2)))**2d0)+2d0/Eur
             df12_hsmdS(3)    = 0d0
             
             df_capmdS(1)     = 2d0*(q_m/(alfa**2d0)+(1d0/3d0)*p)
             df_capmdS(2)     = 2d0*(q_m*(delta-1d0)/(alfa**2d0)+(1d0/3d0)*p)
             df_capmdS(3)     = 2d0*(-q_m*delta/(alfa**2d0)+(1d0/3d0)*p)             

             Call WriVec(ideb, 'df13_hsmdS', df13_hsmdS, 3)                                          
             Call WriVec(ideb, 'df12_hsmdS', df12_hsmdS, 3)                
             Call WriVec(ideb, 'df_capmdS', df_capmdS, 3)         
            
             !Cálculo df13_hs/dq, df12_hs/dq y df_cap/dq en Prs (q1=gamma_p  -para f13_hs y f12_hs-, q2=p_p)
             df13_hsmdq      = -1d0
             df12_hsmdq      = -1d0             
             df_capmdq       = -2d0*p_p 

             Call WriVal(ideb, 'df13_hsmdq', df13_hsmdq)                    
             Call WriVal(ideb, 'df12_hsmdq', df12_hsmdq)                    
             Call WriVal(ideb, 'df_capmdq', df_capmdq)                    
             
             !Cálculo h en Prs (dq=dlam1*h1 -para f13_hs y f12_hs- y dq2=dlam2*h2)
             xh_hs            = 1d0
             xh_cap           = 2d0*xH*p*(max(Prs(3)+c*cos(phi)/sin(phi), c*cos(phi)/sin(phi))/(p_ref+c*cos(phi)/sin(phi)))**xm

             Call WriVal(ideb, 'xh_hs', xh_hs)                    
             Call WriVal(ideb, 'xh_cap', xh_cap)                    
             
             !Cálculo dg13_hs/dS, dg12_hs/dS y dg_cap/dS en Prs
             sin_phi_m        = (Prs(1)-Prs(3))/(Prs(1)+Prs(3)-2d0*c*cos(phi)/sin(phi))
             sin_phi_cv       = (sin(phi)-sin(psi))/(1d0-(sin(phi))*(sin(psi)))
             If (sin_phi_m .LT. (3d0/4d0)*sin(phi)) Then
                 sin_psi_m = 0d0
             ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .GT. 0d0) Then
                 sin_psi_m = max((sin_phi_m-sin_phi_cv)/(1d0-sin_phi_m*sin_phi_cv),0d0)
             ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .LE. 0d0) Then
                 psi_m = psi
                 sin_psi_m = sin(psi_m)
             ElseIf (phi .EQ. 0d0) Then
                 sin_psi_m = 0d0
             End If
             
             dg13_hsmdS(1)   = (1d0/2d0)*(1d0-sin_psi_m)  
             dg13_hsmdS(2)   = 0d0
             dg13_hsmdS(3)   = -(1d0/2d0)*(1d0+sin_psi_m)         

             dg12_hsmdS(1)   = (1d0/2d0)*(1d0-sin_psi_m)  
             dg12_hsmdS(2)   = -(1d0/2d0)*(1d0+sin_psi_m)
             dg12_hsmdS(3)   = 0d0               
              
             dg_capmdS(1)   = df_capmdS(1)     
             dg_capmdS(2)   = df_capmdS(2)     
             dg_capmdS(3)   = df_capmdS(3)     

             Call WriVec(ideb, 'dg13_hsmdS', dg13_hsmdS, 3)                                          
             Call WriVec(ideb, 'dg12_hsmdS', dg12_hsmdS, 3)                                      
             Call WriVec(ideb, 'dg_capmdS', dg_capmdS, 3)                                         
             
             !Cálculo d2g13_hs/dS2, d2g12_hs/dS2 y d2g_cap/dS2 en Prs
             Call MZeroR(d2g13_hsmdS2, 9) 
             Call MZeroR(d2g12_hsmdS2, 9)              
            
             d2g_capmdS2(1,1) = 2d0*(1d0/(alfa**2d0)+1d0/9d0)
             d2g_capmdS2(1,2) = 2d0*((delta-1d0)/(alfa**2d0)+1d0/9d0)
             d2g_capmdS2(1,3) = -2d0*(delta/(alfa**2d0)-1d0/9d0)
             d2g_capmdS2(2,1) = d2g_capmdS2(1,2)    
             d2g_capmdS2(2,2) = 2d0*(((delta-1d0)**2d0)/(alfa**2d0)+1d0/9d0)   
             d2g_capmdS2(2,3) = -2d0*(((delta-1d0)*delta)/(alfa**2d0)-1d0/9d0)   
             d2g_capmdS2(3,1) = d2g_capmdS2(1,3)   
             d2g_capmdS2(3,2) = d2g_capmdS2(2,3)   
             d2g_capmdS2(3,3) = 2d0*((delta**2d0)/(alfa**2d0)+1d0/9d0)   
             
             !Cálculo d2g13_hs/dSdq, d2g12_hs/dSdq y d2g_cap/dSdq en Prs
             Call MZeroR(d2g13_hsmdSdq, 3) 
             Call MZeroR(d2g12_hsmdSdq, 3) 
             Call MZeroR(d2g_capmdSdq, 3)              
             
             !Cálculo del multiplicador plástico 
             Call MatMat(D, 3, d2g13_hsmdS2, 3, 3, 3, 3, D_dg13_2, 3)
             Call MatMat(D, 3, d2g12_hsmdS2, 3, 3, 3, 3, D_dg12_2, 3)             
             Call MatMat(D, 3, d2g_capmdS2, 3, 3, 3, 3, D_dgcap_2, 3)                          
             HH_inv = Unity + xLam13_hs*D_dg13_2 + xLam12_hs*D_dg12_2 + xLam_cap*D_dgcap_2
             Call MatInvPiv(HH_inv, HH, 3)
             Call AddVec(dg13_hsmdS, d2g13_hsmdSdq, 1d0, xLam13_hs*xh_hs, 3, theta13) 
             Call AddVec(dg12_hsmdS, d2g12_hsmdSdq, 1d0, xLam12_hs*xh_hs, 3, theta12) 
             Call AddVec(dg_capmdS, d2g_capmdSdq, 1d0, xLam_cap*xh_cap, 3, thetacap)                      
             Call MatMat(HH, 3, df13_hsmdS, 3, 3, 1, 3, HH_df13,3)
             Call MatMat(HH, 3, df12_hsmdS, 3, 3, 1, 3, HH_df12,3)             
             Call MatMat(HH, 3, df_capmdS, 3, 3, 1, 3, HH_dfcap,3)             
             Call MatMat(D, 3, theta13, 3, 3, 1, 3, D_theta13,3)
             Call MatMat(D, 3, theta12, 3, 3, 1, 3, D_theta12,3)             
             Call MatMat(D, 3, thetacap, 3, 3, 1, 3, D_thetacap,3)              
             Call DInProd(D_theta13, HH_df13, 3, D_theta13_HH_df13) 
             Call DInProd(D_theta12, HH_df13, 3, D_theta12_HH_df13)
             Call DInProd(D_thetacap, HH_df13, 3, D_thetacap_HH_df13)
             Call DInProd(D_theta13, HH_df12, 3, D_theta13_HH_df12)                           
             Call DInProd(D_theta12, HH_df12, 3, D_theta12_HH_df12)    
             Call DInProd(D_thetacap, HH_df12, 3, D_thetacap_HH_df12)    
             Call DInProd(D_theta13, HH_dfcap, 3, D_theta13_HH_dfcap)    
             Call DInProd(D_theta12, HH_dfcap, 3, D_theta12_HH_dfcap)    
             Call DInProd(D_thetacap, HH_dfcap, 3, D_thetacap_HH_dfcap)    
             xOmega11  = D_theta13_HH_df13-df13_hsmdq*xh_hs                     
             xOmega12  = D_theta12_HH_df13
             xOmega13  = D_thetacap_HH_df13
             xOmega21  = D_theta13_HH_df12                              
             xOmega22  = D_theta12_HH_df12-df12_hsmdq*xh_hs
             xOmega23  = D_thetacap_HH_df12                              
             xOmega31  = D_theta13_HH_dfcap                              
             xOmega32  = D_theta12_HH_dfcap                              
             xOmega33  = D_thetacap_HH_dfcap-df_capmdq*xh_cap                              
             Call MatMat(D, 3, dg13_hsmdS, 3, 3, 1, 3, D_dg13, 3)       
             Call MatMat(D, 3, dg12_hsmdS, 3, 3, 1, 3, D_dg12, 3)                      
             Call MatMat(D, 3, dg_capmdS, 3, 3, 1, 3, D_dgcap, 3)               
             Call AddVec(D_dg13, D_dg12, -xLam13_hs, -xLam12_hs, 3, Prs_res1)
             Call AddVec(Prs_res1, D_dgcap, 1d0, -xLam_cap, 3, Prs_res2)       
             Call AddVec(Prs_E, Prs_res2, 1d0, 1d0, 3, Prs_res3)
             Call AddVec(Prs, Prs_res3, 1d0,-1d0, 3, r)         
             Call DInProd(HH_df13, r, 3, r_HH_df13)                           
             Call DInProd(HH_df12, r, 3, r_HH_df12)             
             Call DInProd(HH_dfcap, r, 3, r_HH_dfcap)                  
             xOmegaf13  = f13_hs - r_HH_df13
             xOmegaf12  = f12_hs - r_HH_df12             
             xOmegafcap = f_cap - r_HH_dfcap         
             
             !Resolución sistema y calculo dxLam13_hs y dxLam12_hs            
             xMat(1,1) = xOmega11
             xMat(1,2) = xOmega12
             xMat(1,3) = xOmega13         
             xMat(2,1) = xOmega21         
             xMat(2,2) = xOmega22         
             xMat(2,3) = xOmega23         
             xMat(3,1) = xOmega31         
             xMat(3,2) = xOmega32         
             xMat(3,3) = xOmega33         
             f_b(1) = xOmegaf13
             f_b(2) = xOmegaf12         
             f_b(3) = xOmegafcap         
         
             Call MatInvPiv(xMat, xMat_inv, 3)
             Call MatMat(xMat_inv, 3, f_b, 3, 3, 1, 3, d_lam, 3) 
          
             dxLam13_hs = d_lam(1)
             dxLam12_hs = d_lam(2)
             dxLam_cap  = d_lam(3)
             
             xLam13_hs = xLam13_hs + dxLam13_hs         
             xLam12_hs = xLam12_hs + dxLam12_hs         
             xLam_cap  = xLam_cap + dxLam_cap                               
             
             Call WriVal(ideb, 'dxLam13_hs', dxLam13_hs)                                 
             Call WriVal(ideb, 'dxLam12_hs', dxLam12_hs)                                          
             Call WriVal(ideb, 'dxLam_cap', dxLam_cap)                                              
             Call WriVal(ideb, 'xLam13_hs', xLam13_hs)                                 
             Call WriVal(ideb, 'xLam12_hs', xLam12_hs)                                          
             Call WriVal(ideb, 'xLam_cap', xLam_cap)                   
             
             !Actualizo variables
             Call AddVec(D_dg13, D_dg12, -xLam13_hs, -xLam12_hs, 3, Prs_2sup) 
             Call AddVec(Prs_2sup, D_dgcap, 1d0, -xLam_cap, 3, Prs_3sup)          
             Call AddVec(Prs_E, Prs_3sup, 1d0, 1d0, 3, Prs)
             gamma_p = gamma_p_E + xLam13_hs*xh_hs
             p_p = p_p_E + xLam_cap*xh_cap

             Call WriVal(ideb, 'gamma_p', gamma_p)                                 
             Call WriVal(ideb, 'p_p', p_p)                                              
             
             !Calculo el valor de f13_hs, f12_hs y f_cap en (k+1) en Prs
             qf1         = max(c*cos(phi)/sin(phi)+Prs(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))!OJO, Prs(3) debe ser mayor que c*cos(phi)/sin(phi) para que E sean reales
             qf2         = (2d0*sin(phi)/(1d0-sin(phi)))
             qf          = qf1*qf2
             qa          = qf/R_f
             E_fac1      = max(c*cos(phi)+Prs(3)*sin(phi), c*cos(phi))
             E_fac2      = (c*cos(phi)+p_ref*sin(phi))
             E50         = E50_ref*(E_fac1/E_fac2)**xm
             Ei          = 2d0*E50/(2d0-R_f)
             Eur         = Eur_ref*(E_fac1/E_fac2)**xm          

             f13_hs_fac1 = ((Prs(1)-Prs(3))/(qa-(Prs(1)-Prs(3))))
             f13_hs_fac2 = 2d0*(Prs(1)-Prs(3))/Eur
             f13_hs      = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p

             f12_hs_fac1 = ((Prs(1)-Prs(2))/(qa-(Prs(1)-Prs(2))))
             f12_hs_fac2 = 2d0*(Prs(1)-Prs(2))/Eur
             f12_hs      = (2d0*qa/Ei)*f12_hs_fac1-f12_hs_fac2-gamma_p
                      
             q_m         = Prs(1)+(delta-1d0)*Prs(2)-delta*Prs(3)
             p           = (Prs(1)+Prs(2)+Prs(3))/3d0          
         
             f_cap       = (q_m**2d0)/(alfa**2d0)+p**2d0-p_p**2d0                     
             
             Call WriIvl(ideb, 'it', it)                     
             
             !Actualización iteración
             it = it+1
             
         End do
            
         !Actualización variables de estado
         StVar(1) = gamma_p
         StVar(2) = min((-1d0)*p_p, (-1d0)*p_p_E)

         Call WriVal(ideb, 'xLam13_hs', xLam13_hs)                
         Call WriVal(ideb, 'xLam12_hs', xLam12_hs)    
         Call WriVal(ideb, 'xLam_cap', xLam_cap)             
                  
         If (xLam13_hs .GT. (-1d0)*ttol .AND. &
             xLam12_hs .GT. (-1d0)*ttol .AND. &
             xLam_cap .GT. (-1d0)*ttol .AND. &
             it .LE. itmax) Then
             write(1,*)'Return mapping Ok!'
         Else
             write(1,*)'Numerical problem'
         End If
         
      Case (3) ! 3 : Eje entre f13_hs y f23_hs (Prs(1) = Prs(2) > Prs(3))

         !STARTING POINT
         
         !Cálculo valor de f13_hs, f23_hs y f_cap en Prs_E
         qf1         = max(c*cos(phi)/sin(phi)+Prs_E(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))!OJO, Prs(3) debe ser mayor que c*cos(phi)/sin(phi) para que E sean reales
         qf2         = (2d0*sin(phi)/(1d0-sin(phi)))
         qf          = qf1*qf2
         qa          = qf/R_f
         E_fac1      = max(c*cos(phi)+Prs_E(3)*sin(phi), c*cos(phi))
         E_fac2      = (c*cos(phi)+p_ref*sin(phi))
         E50         = E50_ref*(E_fac1/E_fac2)**xm
         Ei          = 2d0*E50/(2d0-R_f)
         Eur         = Eur_ref*(E_fac1/E_fac2)**xm          
         
         f13_hs_fac1 = ((Prs_E(1)-Prs_E(3))/(qa-(Prs_E(1)-Prs_E(3))))
         f13_hs_fac2 = 2d0*(Prs_E(1)-Prs_E(3))/Eur
         f13_hs      = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p_E

         f23_hs_fac1 = ((Prs_E(2)-Prs_E(3))/(qa-(Prs_E(2)-Prs_E(3))))
         f23_hs_fac2 = 2d0*(Prs_E(2)-Prs_E(3))/Eur
         f23_hs      = (2d0*qa/Ei)*f23_hs_fac1-f23_hs_fac2-gamma_p_E
                  
         xlam_ast    = p_ref/Eoed_ref
         xkap_ast    = 3d0*p_ref*(1-2*xNu_ur)/Eur_ref
         alfa1       = ((1d0-xK0_nc)**2d0)/((1d0+2d0*xK0_nc)**2d0)
         alfa21      = (1d0-xK0_nc)*(1d0-2d0*xNu_ur)*(xlam_ast/xkap_ast-1d0)
         alfa221     = (1d0+2d0*xK0_nc)*(1d0-2d0*xNu_ur)*(xlam_ast/xkap_ast)
         alfa222     = (1d0-xK0_nc)*(1d0+xNu_ur)
         alfa        = 3d0*(alfa1+(alfa21/(alfa221-alfa222)))**0.5d0              
         delta       = (3d0+sin(phi))/(3d0-sin(phi))
         q_m         = Prs_E(1)+(delta-1d0)*Prs_E(2)-delta*Prs_E(3)
         p           = (Prs_E(1)+Prs_E(2)+Prs_E(3))/3d0          

         f_cap       = (q_m**2d0)/(alfa**2d0)+p**2d0-p_p_E**2d0  
         
         !Cálculo tolerancia f13_hs, f23_hs y f_cap en Prs_E
         f13_hs_TOL = ttol*abs(f13_hs)+ttol
         f23_hs_TOL = ttol*abs(f23_hs)+ttol
         f_cap_TOL = ttol*abs(f_cap)+ttol 
         
         !Cálculo df13_hs/dS, df23_hs/dS y df_cap/dS en Prs_E
         df13_hsmdS(1)    = ((2d0*(qa**2d0))/Ei)/((qa-(Prs_E(1)-Prs_E(3)))**2d0)-2d0/Eur
         df13_hsmdS(2)    = 0d0
         df13_hsmdS(3)    = -((2d0*(qa**2d0))/Ei)/((qa-(Prs_E(1)-Prs_E(3)))**2d0)+2d0/Eur 
         
         df23_hsmdS(1)    = 0d0
         df23_hsmdS(2)    = ((2d0*(qa**2d0))/Ei)/((qa-(Prs_E(1)-Prs_E(3)))**2d0)-2d0/Eur
         df23_hsmdS(3)    = -((2d0*(qa**2d0))/Ei)/((qa-(Prs_E(1)-Prs_E(3)))**2d0)+2d0/Eur   
         
         df_capmdS(1)    = 2d0*(q_m/(alfa**2d0)+(1d0/3d0)*p)
         df_capmdS(2)    = 2d0*(q_m*(delta-1d0)/(alfa**2d0)+(1d0/3d0)*p)
         df_capmdS(3)    = 2d0*(-q_m*delta/(alfa**2d0)+(1d0/3d0)*p)         
         
         !Cálculo df13_hs/dq, df23_hs/dq y df_cap/dq en Prs_E (q1=gamma_p  -para f13_hs y f23_hs-, q2=p_p)
         df13_hsmdq      = -1d0
         df23_hsmdq      = -1d0
         df_capmdq       = -2d0*p_p_E            
         
         !Cálculo h en Prs_E (dq=dlam1*h1 -para f13_hs y f23_hs- y dq2=dlam2*h2)
         xh_hs           = 1d0
         xh_cap          = 2d0*xH*p*(max(Prs_E(3)+c*cos(phi)/sin(phi), c*cos(phi)/sin(phi))/(p_ref+c*cos(phi)/sin(phi)))**xm        
         
         !Cálculo dg13_hs/dS, dg23_hs/dS y dg_cap/dS en Prs_E
         sin_phi_m        = (Prs_E(1)-Prs_E(3))/(Prs_E(1)+Prs_E(3)-2d0*c*cos(phi)/sin(phi))
         sin_phi_cv       = (sin(phi)-sin(psi))/(1d0-(sin(phi))*(sin(psi)))
           If (sin_phi_m .LT. (3d0/4d0)*sin(phi)) Then
                sin_psi_m = 0d0
           ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .GT. 0d0) Then
                sin_psi_m = max((sin_phi_m-sin_phi_cv)/(1d0-sin_phi_m*sin_phi_cv),0d0)
           ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .LE. 0d0) Then
                psi_m = psi
                sin_psi_m = sin(psi_m)
           ElseIf (phi .EQ. 0d0) Then
                sin_psi_m = 0d0
           End If
           
         dg13_hsmdS(1)   = (1d0/2d0)*(1d0-sin_psi_m)  
         dg13_hsmdS(2)   = 0d0
         dg13_hsmdS(3)   = -(1d0/2d0)*(1d0+sin_psi_m)         

         dg23_hsmdS(1)   = 0d0
         dg23_hsmdS(2)   = (1d0/2d0)*(1d0-sin_psi_m)  
         dg23_hsmdS(3)   = -(1d0/2d0)*(1d0+sin_psi_m)          

         dg_capmdS(1)   = df_capmdS(1) 
         dg_capmdS(2)   = df_capmdS(2) 
         dg_capmdS(3)   = df_capmdS(3) 
         
         !Cálculo d2g13_hs/dS2, d2g23_hs/dS2 y d2g_cap/dS2 en Prs_E
         Call MZeroR(d2g13_hsmdS2, 9) 
         Call MZeroR(d2g23_hsmdS2, 9)          
         
         d2g_capmdS2(1,1) = 2d0*(1d0/(alfa**2d0)+1d0/9d0)
         d2g_capmdS2(1,2) = 2d0*((delta-1d0)/(alfa**2d0)+1d0/9d0)
         d2g_capmdS2(1,3) = -2d0*(delta/(alfa**2d0)-1d0/9d0)
         d2g_capmdS2(2,1) = d2g_capmdS2(1,2)    
         d2g_capmdS2(2,2) = 2d0*(((delta-1d0)**2d0)/(alfa**2d0)+1d0/9d0)   
         d2g_capmdS2(2,3) = -2d0*(((delta-1d0)*delta)/(alfa**2d0)-1d0/9d0)   
         d2g_capmdS2(3,1) = d2g_capmdS2(1,3)   
         d2g_capmdS2(3,2) = d2g_capmdS2(2,3)   
         d2g_capmdS2(3,3) = 2d0*((delta**2d0)/(alfa**2d0)+1d0/9d0)   
         
         !Cálculo d2g13_hs/dSdq, d2g23_hs/dSdq y d2g_cap/dSdq en Prs_E
         Call MZeroR(d2g13_hsmdSdq, 3) 
         Call MZeroR(d2g23_hsmdSdq, 3)  
         Call MZeroR(d2g_capmdSdq, 3)          
         
         !Cálculo de los multiplicadores plásticos 
         Call MatMat(D, 3, dg13_hsmdS, 3, 3, 1, 3, D_dg13, 3) 
         Call DInProd(df13_hsmdS, D_dg13,3, df13_D_dg13)
         xOmega11 = df13_D_dg13-df13_hsmdq*xh_hs        
         Call MatMat(D,3,dg23_hsmdS,3,3,1,3,D_dg23,3) 
         Call DInProd(df13_hsmdS, D_dg23,3, df13_D_dg23)
         xOmega12 = df13_D_dg23     
         Call MatMat(D,3,dg_capmdS,3,3,1,3,D_dgcap,3) 
         Call DInProd(df13_hsmdS, D_dgcap,3, df13_D_dgcap)
         xOmega13 = df13_D_dgcap              
         Call DInProd(df23_hsmdS, D_dg13,3, df23_D_dg13)
         xOmega21 = df23_D_dg13               
         Call DInProd(df23_hsmdS, D_dg23,3, df23_D_dg23)
         xOmega22 = df23_D_dg23-df23_hsmdq*xh_hs 
         Call DInProd(df23_hsmdS, D_dgcap,3, df23_D_dgcap)
         xOmega23 = df23_D_dgcap
         Call DInProd(df_capmdS, D_dg13,3, dfcap_D_dg13)
         xOmega31 = dfcap_D_dg13
         Call DInProd(df_capmdS, D_dg23,3, dfcap_D_dg23)
         xOmega32 = dfcap_D_dg23
         Call DInProd(df_capmdS, D_dgcap,3, dfcap_D_dgcap)
         xOmega33 = dfcap_D_dgcap-df_capmdq*xh_cap          
                         
         !Resolución sistema y cálculo dxLam13_hs, dxLam23_hs y dxLam_cap
         xMat(1,1) = xOmega11
         xMat(1,2) = xOmega12
         xMat(1,3) = xOmega13         
         xMat(2,1) = xOmega21         
         xMat(2,2) = xOmega22         
         xMat(2,3) = xOmega23         
         xMat(3,1) = xOmega31         
         xMat(3,2) = xOmega32         
         xMat(3,3) = xOmega33         
         f_b(1) = f13_hs
         f_b(2) = f23_hs         
         f_b(3) = f_cap         
         
         Call MatInvPiv(xMat, xMat_inv, 3)
         Call MatMat(xMat_inv, 3, f_b, 3, 3, 1, 3, d_lam, 3) 
         
         dxLam13_hs = d_lam(1)
         dxLam23_hs = d_lam(2)
         dxLam_cap  = d_lam(3)
         
         xLam13_hs = dxLam13_hs         
         xLam23_hs = dxLam23_hs         
         xLam_cap  = dxLam_cap                  
         
         !Actualización variables
         Call AddVec(D_dg13, D_dg23, -xLam13_hs, -xLam23_hs, 3, Prs_2sup) 
         Call AddVec(Prs_2sup, D_dgcap, 1d0, -xLam_cap, 3, Prs_3sup)          
         Call AddVec(Prs_E, Prs_3sup, 1d0, 1d0, 3, Prs)
         gamma_p = gamma_p_E + xLam13_hs*xh_hs
         p_p = p_p_E + xLam_cap*xh_cap
       
         !Cálculo valor de f13_hs, f23_hs y f_cap con la nueva Prs         
         qf1         = max(c*cos(phi)/sin(phi)+Prs(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))!OJO, Prs(3) debe ser mayor que c*cos(phi)/sin(phi) para que E sean reales
         qf2         = (2d0*sin(phi)/(1d0-sin(phi)))
         qf          = qf1*qf2
         qa          = qf/R_f
         E_fac1      = max(c*cos(phi)+Prs(3)*sin(phi), c*cos(phi))
         E_fac2      = (c*cos(phi)+p_ref*sin(phi))
         E50         = E50_ref*(E_fac1/E_fac2)**xm
         Ei          = 2d0*E50/(2d0-R_f)
         Eur         = Eur_ref*(E_fac1/E_fac2)**xm          

         f13_hs_fac1 = ((Prs(1)-Prs(3))/(qa-(Prs(1)-Prs(3))))
         f13_hs_fac2 = 2d0*(Prs(1)-Prs(3))/Eur
         f13_hs      = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p

         f23_hs_fac1 = ((Prs(2)-Prs(3))/(qa-(Prs(2)-Prs(3))))
         f23_hs_fac2 = 2d0*(Prs(2)-Prs(3))/Eur
         f23_hs      = (2d0*qa/Ei)*f23_hs_fac1-f23_hs_fac2-gamma_p
         
         q_m         = Prs(1)+(delta-1d0)*Prs(2)-delta*Prs(3)
         p           = (Prs(1)+Prs(2)+Prs(3))/3d0          
         
         f_cap       = (q_m**2d0)/(alfa**2d0)+p**2d0-p_p**2d0     
         
         !EULER BACKWARD
         
         Do while ((abs(f13_hs) .GT. f13_hs_TOL .OR. &
                   abs(f23_hs) .GT. f23_hs_TOL .OR. &
                   abs(f_cap) .GT. f_cap_TOL) .AND. &
                   it .LE. itmax)
         
             !Construcción matriz D
             !Call Matriz3(Props, Prs, D, IsUndr)                          
             
             !Cálculo df13_hs/dS, df23_hs/dS y df_cap/dS en Prs
             df13_hsmdS(1)    = ((2d0*(qa**2d0))/Ei)/((qa-(Prs(1)-Prs(3)))**2d0)-2d0/Eur
             df13_hsmdS(2)    = 0d0
             df13_hsmdS(3)    = -((2d0*(qa**2d0))/Ei)/((qa-(Prs(1)-Prs(3)))**2d0)+2d0/Eur 
             
             df23_hsmdS(1)    = 0d0
             df23_hsmdS(2)    = ((2d0*(qa**2d0))/Ei)/((qa-(Prs(1)-Prs(2)))**2d0)-2d0/Eur
             df23_hsmdS(3)    = -((2d0*(qa**2d0))/Ei)/((qa-(Prs(1)-Prs(3)))**2d0)+2d0/Eur

             df_capmdS(1)     = 2d0*(q_m/(alfa**2d0)+(1d0/3d0)*p)
             df_capmdS(2)     = 2d0*(q_m*(delta-1d0)/(alfa**2d0)+(1d0/3d0)*p)
             df_capmdS(3)     = 2d0*(-q_m*delta/(alfa**2d0)+(1d0/3d0)*p)             
             
            !Cálculo df13_hs/dq, df23_hs/dq y df_cap/dq en Prs (q1=gamma_p  -para f13_hs y f23_hs-, q2=p_p)
             df13_hsmdq      = -1d0
             df23_hsmdq      = -1d0
             df_capmdq        = -2d0*p_p 
         
            !Cálculo h en Prs (dq=dlam1*h1 -para f13_hs y f23_hs- y dq2=dlam2*h2)
             xh_hs            = 1d0
             xh_cap           = 2d0*xH*p*(max(Prs(3)+c*cos(phi)/sin(phi), c*cos(phi)/sin(phi))/(p_ref+c*cos(phi)/sin(phi)))**xm
             
             !Cálculo dg13_hs/dS, dg23_hs/dS y dg_cap/dS en Prs
             sin_phi_m        = (Prs(1)-Prs(3))/(Prs(1)+Prs(3)-2d0*c*cos(phi)/sin(phi))
             sin_phi_cv       = (sin(phi)-sin(psi))/(1d0-(sin(phi))*(sin(psi)))
             If (sin_phi_m .LT. (3d0/4d0)*sin(phi)) Then
                 sin_psi_m = 0d0
             ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .GT. 0d0) Then
                 sin_psi_m = max((sin_phi_m-sin_phi_cv)/(1d0-sin_phi_m*sin_phi_cv),0d0)
             ElseIf (sin_phi_m .GE. (3d0/4d0)*sin(phi) .AND. psi .LE. 0d0) Then
                 psi_m = psi
                 sin_psi_m = sin(psi_m)
             ElseIf (phi .EQ. 0d0) Then
                 sin_psi_m = 0d0
             End If
             
             dg13_hsmdS(1)   = (1d0/2d0)*(1d0-sin_psi_m)  
             dg13_hsmdS(2)   = 0d0
             dg13_hsmdS(3)   = -(1d0/2d0)*(1d0+sin_psi_m)         

             dg23_hsmdS(1)   = 0d0
             dg23_hsmdS(2)   = (1d0/2d0)*(1d0-sin_psi_m)  
             dg23_hsmdS(3)   = -(1d0/2d0)*(1d0+sin_psi_m)                 
              
             dg_capmdS(1)    = df_capmdS(1)
             dg_capmdS(2)    = df_capmdS(2)
             dg_capmdS(3)    = df_capmdS(3)
                      
             !Cálculo d2g13_hs/dS2, d2g23_hs/dS2 y d2g_cap/dS2 en Prs
             Call MZeroR(d2g13_hsmdS2, 9) 
             Call MZeroR(d2g23_hsmdS2, 9)              
          
             d2g_capmdS2(1,1) = 2d0*(1d0/(alfa**2d0)+1d0/9d0)
             d2g_capmdS2(1,2) = 2d0*((delta-1d0)/(alfa**2d0)+1d0/9d0)
             d2g_capmdS2(1,3) = -2d0*(delta/(alfa**2d0)-1d0/9d0)
             d2g_capmdS2(2,1) = d2g_capmdS2(1,2)    
             d2g_capmdS2(2,2) = 2d0*(((delta-1d0)**2d0)/(alfa**2d0)+1d0/9d0)   
             d2g_capmdS2(2,3) = -2d0*(((delta-1d0)*delta)/(alfa**2d0)-1d0/9d0)   
             d2g_capmdS2(3,1) = d2g_capmdS2(1,3)   
             d2g_capmdS2(3,2) = d2g_capmdS2(2,3)   
             d2g_capmdS2(3,3) = 2d0*((delta**2d0)/(alfa**2d0)+1d0/9d0)   
             
             !Cálculo d2g13_hs/dSdq, d2g23_hs/dSdq y d2g_cap/dSdq en Prs
             Call MZeroR(d2g13_hsmdSdq, 3) 
             Call MZeroR(d2g23_hsmdSdq, 3) 
             Call MZeroR(d2g_capmdSdq, 3)              
             
             !Cálculo del multiplicador plástico 
             Call MatMat(D, 3, d2g13_hsmdS2, 3, 3, 3, 3, D_dg13_2, 3)
             Call MatMat(D, 3, d2g23_hsmdS2, 3, 3, 3, 3, D_dg23_2, 3)             
             Call MatMat(D, 3, d2g_capmdS2, 3, 3, 3, 3, D_dgcap_2, 3)                          
             HH_inv = Unity + xLam13_hs*D_dg13_2 + xLam23_hs*D_dg23_2 + xLam_cap*D_dgcap_2
             Call MatInvPiv(HH_inv, HH, 3)
             Call AddVec(dg13_hsmdS, d2g13_hsmdSdq, 1d0, xLam13_hs*xh_hs, 3, theta13) 
             Call AddVec(dg23_hsmdS, d2g23_hsmdSdq, 1d0, xLam23_hs*xh_hs, 3, theta23) 
             Call AddVec(dg_capmdS, d2g_capmdSdq, 1d0, xLam_cap*xh_cap, 3, thetacap) 
             Call MatMat(HH, 3, df13_hsmdS, 3, 3, 1, 3, HH_df13,3)
             Call MatMat(HH, 3, df23_hsmdS, 3, 3, 1, 3, HH_df23,3)             
             Call MatMat(HH, 3, df_capmdS, 3, 3, 1, 3, HH_dfcap,3)             
             Call MatMat(D, 3, theta13, 3, 3, 1, 3, D_theta13,3)
             Call MatMat(D, 3, theta23, 3, 3, 1, 3, D_theta23,3)             
             Call MatMat(D, 3, thetacap, 3, 3, 1, 3, D_thetacap,3)              
             Call DInProd(D_theta13, HH_df13, 3, D_theta13_HH_df13) 
             Call DInProd(D_theta23, HH_df13, 3, D_theta23_HH_df13)
             Call DInProd(D_thetacap, HH_df13, 3, D_thetacap_HH_df13)
             Call DInProd(D_theta13, HH_df23, 3, D_theta13_HH_df23)                           
             Call DInProd(D_theta23, HH_df23, 3, D_theta23_HH_df23)    
             Call DInProd(D_thetacap, HH_df23, 3, D_thetacap_HH_df23)    
             Call DInProd(D_theta13, HH_dfcap, 3, D_theta13_HH_dfcap)    
             Call DInProd(D_theta23, HH_dfcap, 3, D_theta23_HH_dfcap)    
             Call DInProd(D_thetacap, HH_dfcap, 3, D_thetacap_HH_dfcap)    
             xOmega11  = D_theta13_HH_df13-df13_hsmdq*xh_hs                     
             xOmega12  = D_theta23_HH_df13
             xOmega13  = D_thetacap_HH_df13
             xOmega21  = D_theta13_HH_df23                              
             xOmega22  = D_theta23_HH_df23-df23_hsmdq*xh_hs
             xOmega23  = D_thetacap_HH_df23                              
             xOmega31  = D_theta13_HH_dfcap                              
             xOmega32  = D_theta23_HH_dfcap                              
             xOmega33  = D_thetacap_HH_dfcap-df_capmdq*xh_cap                              
             Call MatMat(D, 3, dg13_hsmdS, 3, 3, 1, 3, D_dg13,3)       
             Call MatMat(D, 3, dg23_hsmdS, 3, 3, 1, 3, D_dg23, 3)                      
             Call MatMat(D, 3, dg_capmdS, 3, 3, 1, 3, D_dgcap, 3)               
             Call AddVec(D_dg13, D_dg23, -xLam13_hs, -xLam23_hs, 3, Prs_res1)
             Call AddVec(Prs_res1, D_dgcap, 1d0, -xLam_cap, 3, Prs_res2)       
             Call AddVec(Prs_E, Prs_res2, 1d0, 1d0, 3, Prs_res3)
             Call AddVec(Prs, Prs_res3, 1d0, -1d0, 3, r)         
             Call DInProd(HH_df13, r, 3, r_HH_df13)                           
             Call DInProd(HH_df23, r, 3, r_HH_df23)             
             Call DInProd(HH_dfcap, r, 3, r_HH_dfcap)                  
             xOmegaf13  = f13_hs - r_HH_df13
             xOmegaf23  = f23_hs - r_HH_df23             
             xOmegafcap = f_cap - r_HH_dfcap         
             
             !Resolución sistema y cálculo dxLam13_hs y dxLam23_hs            
             xMat(1,1) = xOmega11
             xMat(1,2) = xOmega12
             xMat(1,3) = xOmega13         
             xMat(2,1) = xOmega21         
             xMat(2,2) = xOmega22         
             xMat(2,3) = xOmega23         
             xMat(3,1) = xOmega31         
             xMat(3,2) = xOmega32         
             xMat(3,3) = xOmega33         
             f_b(1) = xOmegaf13
             f_b(2) = xOmegaf23         
             f_b(3) = xOmegafcap         
         
             Call MatInvPiv(xMat, xMat_inv, 3)
             Call MatMat(xMat_inv, 3, f_b, 3, 3, 1, 3, d_lam, 3) 
          
             dxLam13_hs = d_lam(1)
             dxLam23_hs = d_lam(2)
             dxLam_cap  = d_lam(3)
             
             xLam13_hs = xLam13_hs + dxLam13_hs         
             xLam23_hs = xLam23_hs + dxLam23_hs         
             xLam_cap  = xLam_cap + dxLam_cap                  
                         
             !Actualización variables
             Call AddVec(D_dg13, D_dg23, -xLam13_hs, -xLam23_hs, 3, Prs_2sup) 
             Call AddVec(Prs_2sup, D_dgcap, 1d0, -xLam_cap, 3, Prs_3sup)          
             Call AddVec(Prs_E, Prs_3sup, 1d0, 1d0, 3, Prs)
             gamma_p = gamma_p_E + xLam13_hs*xh_hs
             p_p = p_p_E + xLam_cap*xh_cap
             
             !Cálculo valor de f13_hs, f23_hs y f_cap en (k+1) en Prs
             qf1         = max(c*cos(phi)/sin(phi)+Prs(3), max(c*cos(phi)/sin(phi)-sigma_trac, 0d0))!OJO, Prs(3) debe ser mayor que c*cos(phi)/sin(phi) para que E sean reales
             qf2         = (2d0*sin(phi)/(1d0-sin(phi)))
             qf          = qf1*qf2
             qa          = qf/R_f
             E_fac1      = max(c*cos(phi)+Prs(3)*sin(phi), c*cos(phi))
             E_fac2      = (c*cos(phi)+p_ref*sin(phi))
             E50         = E50_ref*(E_fac1/E_fac2)**xm
             Ei          = 2d0*E50/(2d0-R_f)
             Eur         = Eur_ref*(E_fac1/E_fac2)**xm          
             
             f13_hs_fac1 = ((Prs(1)-Prs(3))/(qa-(Prs(1)-Prs(3))))
             f13_hs_fac2 = 2d0*(Prs(1)-Prs(3))/Eur
             f13_hs      = (2d0*qa/Ei)*f13_hs_fac1-f13_hs_fac2-gamma_p

             f23_hs_fac1 = ((Prs(2)-Prs(3))/(qa-(Prs(2)-Prs(3))))
             f23_hs_fac2 = 2d0*(Prs(2)-Prs(3))/Eur
             f23_hs      = (2d0*qa/Ei)*f23_hs_fac1-f23_hs_fac2-gamma_p
    
             q_m         = Prs(1)+(delta-1d0)*Prs(2)-delta*Prs(3)
             p           = (Prs(1)+Prs(2)+Prs(3))/3d0          
         
             f_cap       = (q_m**2d0)/(alfa**2d0)+p**2d0-p_p**2d0                     

             Call WriIvl(ideb, 'it', it)                         
             
             !Actualizo iteración
             it = it+1
             
         End do
            
         !Actualizo variables de estado
         StVar(1) = gamma_p
         StVar(2) = min((-1d0)*p_p, (-1d0)*p_p_E)
 
        Call WriVal(ideb, 'xLam13_hs', xLam13_hs)                
        Call WriVal(ideb, 'xLam23_hs', xLam23_hs)    
        Call WriVal(ideb, 'xLam_cap', xLam_cap)             
         
         If (xLam13_hs .GT. (-1d0)*ttol .AND. &
             xLam23_hs .GT. (-1d0)*ttol .AND. &
             xLam_cap .GT. (-1d0)*ttol .AND. &
             it .LE. itmax) Then
             write(1,*)'Return mapping Ok!'
         Else
             write(1,*)'Numerical problem'
         End If


      Case Default
             write(1,*) ' incorrect iArea in HS_plast'
             Stop 
      End Select
      
      
      End ! HS_CAP_plast
____________________________________________________________________________________________________
 
 
CAP_PLAST.F90 , HS_PLAST.F90, HS_TRAC_PLAST.F90, MATRIZ6.f90,  ETC... (OTHER SUBROUTINES)
_____________________________________________________________________________________________________
 
ETC...
ETC...
ETC...
_____________________________________________________________________________________________________
 
 
0 Kudos
mecej4
Honored Contributor III
1,149 Views

OK, using the link that you gave in the second paragraph of #5, I downloaded the Zip file http://kb.plaxis.nl/sites/kb.plaxis.nl/files/kb-downloads/udsm_examples_2012.zip. I built the 32-bit DLL using Intel Fortran 16 (32-bit) after fixing a couple of errors in the batch and options files that were in the Zip:

  • In file ifcomp.bat, surround the replace > by ">", i.e., surround the > character with double quotes.
  • In file comp_rel, change /optimize:4 to /O2

What they have done is to have the batch file first put the compiler-specific directive/declarations for the DLL entry point user_mod into a one-line file (impexp)using the echo command in the batch file. This line is then included in the source file usrmod.for. You do not have this INCLUDE statement in your source file in #5. Without it, the exported name will not be correctly seen when Plaxis attempts to load the DLL.

Everything that I have seen about Plaxis today tells me that they expect users to use one compiler to build the DLL. All that they care about is that the DLL routine be callable by Delphi calling protocol. Therefore, you can use the batch file ifcomp.bat to build the DLL (assuming that your copy of Plaxis wants to call a 32-bit DLL), and ignore the other batch files.

0 Kudos
Javier_C_3
Beginner
1,149 Views

As I'm using Plaxis in a 64 bit environtment, it's necessary to compile the corresponding 64-bit .dll file too (using ifcomp64.bat).

I've done what you pointed in the previous post, but Plaxis still not worrking (I mean Plaxis can't read all necessary information from the .dll's, because it doesn't let me to introduce model parameters in its input page).

In 64 bit dll compilation appeared next message: 
">copy impivf impexp  
System cannot find the file specified"
(however, this message has appeared always since I started compiling weeks ago) 

Yes, I erased "include 'impexp'" statement, as you can see in my Usemod.f90 (copy in #5). What I do to compile using Intel is to activate the line:  !DEC$ ATTRIBUTES DLLExport, StdCall, Reference :: User_Mod  in Usermod.f90.

I think that until now, there were some necessary information for Plaxis missed in the dynamic libraries example_ivf.dll and example_ivf64.dll (due an error in compilation? maybe the error I pointed out before in 64bit dll compilation?).
And  maybe Plaxis found this information in the dynamic library example_lah.dll (compiled with LF95).
I guess this because once add example_lah.dll to Plaxis folder, it work well (I mean that let me introduce the model parameters; this is the only thing that change adding example_lah.dll, but i'ts primordial). T
hat's why I initially thoght that both compilers were necessary.

This supposed missed information in example_ivf.dll and example_ivf64.dll due to a possible error in compilation, might be related with user_add.f90 / useradddf.f90 , which are the responsible of parameter input in Plaxis. I've copy user_add.f90 and useradddf.f90 code (and others .f90) in my last post.

Now, I have the following orders in make files, but still not working:

IFCOMP.BAT
echo !DEC$ ATTRIBUTES DLLExport, StdCall, reference :: User_Mod">"impexp
ifort @comp_rel  UsrMod.f90 UsrAddDF.f90 usrlib.f90 okmsgivf.f90 /dll /exe:example_ivf

IMPEXP
: User_Mod

COMP_REL
/real_size:64
/fpconstant
/O2
/static
/assume:noaccuracy_sensitive
/assume:byterecl
/assume:buffered_io
/notraceback
/debug:none
/check:none
/heap-arrays:1
/fpp
/recursive

IFCOMP64.BAT
REM compile 64-bit DLL: make sure you're in a 64-bit development environment
copy impivf impexp
ifort /static  UsrMod.f90 usrlib.f90 okmsgivf.f90 /dll /exe:example_ivf64

 

 

And, as I said before, there must be some error/missed orders in the previous files (IFCOMP.BAT, IMPEXP, COMP_REL, IFCOMP64.BAT) that is contained in the next one:

LFCOMP.BAT
echo       DLL_Export User_Mod
lf95 -c usrlib.f90
lib /out:lfusrlib.lib usrlib.obj
lf95 -c -ml winapi okmsglah.f90
lf95 -win -dll  UsrMod.f90 Usr_Add.f90 okmsglah.obj -lib lfusrlib.lib -ml bd -out:example_lah.dll

 

Thank you so much for your support.

 

0 Kudos
mecej4
Honored Contributor III
1,149 Views

Page 175 of http://www.plaxis.nl/files/files/2D-3-Material-Models.pdf says:

In order to compile as 64-bit, you need both a 32-bit compiled 'USRMOD.DLL' and a
64-bit compiled ' USRMOD64.DLL' file in the usdm folder under the PLAXIS program
directory. The last one only needs to contain the 'User_Mod' subroutine.

I take this to imply that you should produce two versions of your DLL, one 32-bit and the other 64-bit. For the example code provided in http://kb.plaxis.nl/sites/kb.plaxis.nl/files/kb-downloads/udsm_examples_2012.zip, all that I had to do was

Open a 32-bit IFort command window, create a new directory and extract the examples zip file and correct comp_rel  as in #6. Delete the extracted DLLs or put them into a different directory. Issue the commands

echo !DEC$ ATTRIBUTES DLLExport, StdCall, reference :: User_Mod > impexp
ifort @comp_rel  UsrMod.for UsrAddDF.for usrlib.for okmsgivf.for /dll /exe:example_ivf

Open a 64-bit IFort command window, change to the same working directory as in the previous step, and issue the commands

echo !DEC$ ATTRIBUTES DLLExport :: User_Mod > impexp
ifort @comp_rel  UsrMod.for usrlib.for okmsgivf.for /dll /exe:example_ivf64

Note that I did not use the troublesome batch files with their redirected echo commands, nor did I use the one-line file impivf that is in the Plaxis examples zip but you do not have (or have erased).

The two DLLs were produced. They are not the same as the DLLs in the Plaxis examples zip, probably because they used an older version of the compiler.

To avoid confusion, I suggest that you re-do what I have described in #8 in a separate directory into which you have extracted the Zip file from Plaxis. Once this works on your system, you can move to the directory with your own sources and repeat the procedure but with your sources instead.

0 Kudos
Javier_C_3
Beginner
1,149 Views

I did what you post on #8 with Plaxis Zip (erasing the line with " include 'impexp' " in usrmod.for, changing optimize:4 to O2 in comp_rel, and issuing the comands you've written in #8 in my IA-32 and IA-64 command windows

Dll's are created, but Plaxis still not working, because it doesn't let me introduce model parameters once opened Plaxis Console and selected "Used-defined model" in Material Model type.

As it's explained in the manual (pp. 172-174), Plaxis uses the procedures in UsrAddDF.for to activate user-interface functions (this is need for model parameters input and its recognition). I think that something is going wrong related with this. 

Each procedure has it's own directive (I guess all is correct ¿?):

!DEC$ ATTRIBUTES DLLExport, StdCall, reference ::  GetModelCount
!DEC$ ATTRIBUTES DLLExport, StdCall, reference ::  GetModelName
!DEC$ ATTRIBUTES DLLExport, StdCall, reference ::  GetParamCount
!DEC$ ATTRIBUTES DLLExport, StdCall, reference ::  GetParamName
!DEC$ ATTRIBUTES DLLExport, StdCall, reference ::  GetParamUnit
!DEC$ ATTRIBUTES dllexport, stdcall, reference :: GetStateVarCount
!DEC$ ATTRIBUTES dllexport, stdcall, reference :: GetStateVarName
!DEC$ ATTRIBUTES dllexport, stdcall, reference :: GetStateVarUnit


I've tried also some variations of what you've told me to do:

- Activate the following line in UsrMod.for  !DEC$ ATTRIBUTES DLLExport, StdCall, Reference :: User_Mod . Try compilation.
- Activate the following line in UsrMod.for: include 'usr_add.for' ! give parameter names and units and erase UsrAddDF.for from issue commands. Try compilation.


But none of them fix the problem.

0 Kudos
mecej4
Honored Contributor III
1,149 Views

Javier C. wrote:

I did what you post on #8 with Plaxis Zip (erasing the line with " include 'impexp' " in usrmod.for, changing optimize:4 to O2 in comp_rel, and issuing the comands you've written in #8 in my IA-32 and IA-64 command windows) 

You should not do that. The include 'impexp' line is essential for making the Plaxis-DLL interface work. What gets included (just one line) has to be different for 32-bit and 64-bit, because the calling conventions are different. 

If you prefer, you could make two copies of usrmod.for, say usrmod32.for and usrmod64.for, in which you replace the include 'impexp' line by 

!DEC$ ATTRIBUTES DLLExport, StdCall, reference :: User_Mod 

in usrmod32.for and by

!DEC$ ATTRIBUTES DLLExport :: User_Mod

in usrmod64.for. When you build the DLLs, use usrmod32.for for the 32-bit DLL, and usrmod64.for for the 64-bit DLL.

Beyond that, I cannot comment about the !DEC$ lines that you showed in #9 for the other routines because I do not have Plaxis and I do not know what its requirements are. The documentation seems to indicate that those routines are not needed in the 64-bit DLL, when it says "The last one (USRMOD64.DLL) only needs to contain the 'User_Mod' subroutine.".

0 Kudos
Javier_C_3
Beginner
1,149 Views

Oh, ok, I'm Sorry.

Now I have two copies of UsrMod (UsrMod32.for and UsrMod64.for) with its corresponding lines (post #10) containing !DEC$ directives in substitution of include 'impexp' (I understand now which was the function of line include 'ìmpexp' UsrMod.for, thank you).

Then, I write issues commands:

ifort @comp_rel  UsrMod32.for UsrAddDF.for usrlib.for okmsgivf.for /dll /exe:example_ivf

for UsrMod32.for for 32-bit DLL and

ifort @comp_rel  UsrMod64.for usrlib.for okmsgivf.for /dll /exe:example_ivf64

for UsrMod64.for for 64-bit DLL

Then two .dll are created (example_ivf.dll and example_ivf64.dll). I have copy them in Plaxis Folder.

After that, I open Plaxis, but when I'm going to introduce Model Input Parameters, Plaxis doesn't shown anything.

If all was correct, in the table of Input Parameters of the Model (while executing Plaxis) it should appear the name of the .dll that contains this information, the name of the model, the name of the parameters, their units, etc. All this information is provided by UsrAddDF.for. 

UsrAddDF.for contains the code and directives that let Plaxis "read" how many parameters have the model/s that we are defining, what are their names and units, etc. (I don't know much more about UsrAddDF.for; in my code, I just adapted UsrAddDF from the one in Plaxis ZIP). 

How I said before, I solved this using LF95 (now I can't, since License has expired). To compile with LF95, I only changed nex line

 !DEC$ ATTRIBUTES DLLExport, StdCall, reference :: User_Mod 

in UsrMod by

 DLL_Export User_Mod

Then, after execute lfcomp.bat with LF95, the file example_lah.dll was created and once added together with the files example_ivf.dll and example_ivf64.dll (compiled with Intel) in Plaxis Folder, all works well (Plaxis could read parameters, etc.)

Maybe some clue is in lfcomp.bat (used in LF95 compilation).

Thank you for all your time and comments, and sorry for my ignorance in this issue.

0 Kudos
Javier_C_3
Beginner
1,149 Views

I think I was wrong in one thing. What made Plaxis run until now when using the option of User Soil Defined Model was example_lah.dll compiled with LF95. I've tried to put and old version of this dll alone in Plaxis Folder (without example_ivf.dll and example_ivf64.dll), and it works well.

In any case, I'm sure that two dll compiled with Intel will finally works. 

0 Kudos
Reply