How to implement new soil models in CRISP using the Material Models Interface

The Material Models Interface (MMI) was originally developed by Dr Andrew Chan (Birmingham University) to help add new soil models into CRISP more easily.

There are four important stages in CRISP which call the MMI. These are

A. Reading in the material properties (subroutine MSUB1)

B. Initialising in-situ stresses, yield surface and other initial parameters (subroutine MSINSIT)

C. Assembling the stiffness (subroutine FRONTZ for symmetric stiffness or FRONTU for unsymmetric stiffness)

D. Evaluating the stresses and updating various parameters (subroutine UPOUT)

Stages C and D are passed through within an incremental/iterative solution, while stages A and B are called only at the start of the program run.

The MMI is represented by subroutine CRSM2D and is called through each of the stages above with a flag which indicates at which stage the program is.

The flags are as follows:

ISWDP=1 Program is reading properties (stage A above)

ISWDP=5 Program is initialising stresses (stage B above)

ISWDP=3 or 4 Program is assembling stiffness (stage C above).

ISWDP=2 Program is evaluating stresses (stage D above)

ISWD=3 is for symmetric stiffness, ISWD=4 is for unsymmetric stiffness which would require the special Frontal solver in CRISP for unsymmetric stiffness

To implement a new model, you would need to write your code using an IF condition which uses the switches above. You would need to chose a model ID for your new model. A list of existing model IDs is shown below. You would need to avoid using existing model IDs.

For example you would have the following statement added in subroutine CRSM2D:

if(model.eq.159) call mymodel(propd,lprpd,ISWDP,dstre,dstan,dmatx,gausm


Then create a separate file which includes your model as follows

Subroutine mymodel(propd,lprpd,iswdp,dstre,dstan,dmatx,ester




!.... input/output variables:

dimension propd(lprpd),estre(nstre),param(lpara)

!.... output variables:

dimension dstre(nstre),dmatx(nstre,nstre)

!.some more dimension statements here for local arrays



if(ISWD.eq.1) then


! Read your material properties here. These would be read from record D in MPD file


elseif(ISWD.eq.5) then


.initialise your in-situ stress parameters (such as size of yield stress, etc)


elseif(ISWD.eq.3) then



assemble your stiffness (ie D matrix) here



elseif(ISWD.eq.4) then


.evaluate your stresses, check for yield, etc




The two common blocks above are important as they contain the channel number for reading from the MPD where you would have your material properties (in Record D of the MPD file) as well as other information needed for reading input.

The array ester contains the stresses Sigma_x, Sigma_y, Sigma_z, Tau_xy, Tau_yz, Tau_zx

The array param contains the pore water pressure in the first component (ie param(1)). Subsequent component contain other stress related parameters such as Yield locus and voids ratio in the case of cam-clay models. The arrays ester and param cover all the integration Gauss points of all the elements. Both of these arrays are equivalent to the old VARINT array in CRISP

Current models in CRISP are:

Model ID Model Name
101 Elastic Isotropic
102 Elastic, linear variation with depth
103 Elastic Anisotropic
121 Duncan and Chang Hyperbolic model
122 Spare (for Lade model)
123 Gunn small strain model for undrained problems
124 New small strain model for consolidation analysis
141 von Mises
142 Tresca
143 Drucker Prager
144 Mohr Coulomb
145 Mohr Coulomb elastic-rigid plastic with non-associated flow
155 Advanced Mohr Coulomb hardening model with non-associated flow
157 spare
161 Original Cam Clay
162 Modified Cam Clay
163 Schofield Cam Clay
164 Three Surface Kinematic Hardening model (3-SKH)
165 New 3-SKH
166 MIT-E3 (incomplete)
181 Structural model - Bar
182 Structural model - Beam
191 Structural model - Slip

If you require further assitance please contact Amir