1C --------------------------------------------------------------------
2C MISSING Simulates missing observations
3C Developed by A.Rossi, C.Planas and G.Fiorentini
4C
5C Copyright (C) 2010-2014 European Commission
6C
7C This file is part of Program DMM
8C
9C DMM is free software developed at the Joint Research Centre of the
10C European Commission: you can redistribute it and/or modify it under
11C the terms of the GNU General Public License as published by
12C the Free Software Foundation, either version 3 of the License, or
13C (at your option) any later version.
14C
15C DMM is distributed in the hope that it will be useful,
16C but WITHOUT ANY WARRANTY; without even the implied warranty of
17C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18C GNU General Public License for more details.
19C
20C You should have received a copy of the GNU General Public License
21C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
22C --------------------------------------------------------------------
23	SUBROUTINE MISSING(yk,ny,nz,nx,nu,ns,nt,nmis,theta,S,STATE,pdll,ykmis)
24
25	USE dfwin
26	INTERFACE
27	 SUBROUTINE DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
28	 INTEGER ny,nz,nx,nu,ns(6),nt
29	 DOUBLE PRECISION theta(nt)
30	 DOUBLE PRECISION c(ny,max(1,nz),ns(1)),H(ny,nx,ns(2)),
31	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6))
32	 END SUBROUTINE
33	END INTERFACE
34	CHARACTER*1 fittizia
35	POINTER (pdll,fittizia)
36	POINTER (pdesign,DESIGN)
37
38C INPUT
39	INTEGER ny,nz,nx,nu,nt,ns(6),nmis
40	DOUBLE PRECISION yk(ny+nz),theta(nt),STATE(nx)
41
42C OUTPUT
43	DOUBLE PRECISION ykmis(nmis)
44C LOCALS
45      INTEGER S(6),I,J,K,IFAIL,NIYK(ny)
46	DOUBLE PRECISION U(nu)
47      DOUBLE PRECISION gennor
48	DOUBLE PRECISION,ALLOCATABLE::R(:,:,:),c(:,:,:),H(:,:,:),
49	1 G(:,:,:),a(:,:),F(:,:,:)
50
51	ALLOCATE(R(nx,nu,ns(6)),c(ny,max(nz,1),ns(1)),H(ny,nx,ns(2)),
52	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)))
53	pdesign = getprocaddress(pdll, "design_"C)
54	CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
55
56C NIYK = not(IYK)
57	K = 0
58	DO 10 J = 1,ny
59	 IF(yk(J).EQ.-99999.D0) THEN
60        K = K+1
61	  NIYK(K) = J
6210	 ENDIF
63
64C SAMPLING U
65	IFAIL = -1
66	U(1:nu) = 0.D0
67	DO 20 I = 1,nu
68c	CALL G05EAF(U(I),1,1.D0,1,1.D-14,WORKU,3,IFAIL)
69c20   CALL G05EZF(U(I),1,WORKU,3,IFAIL)
7020    U(I) = gennor(0.D0,1.D0)
71
72
73C DRAW yk ~ f(yk|x,S,zk,theta)
74	DO 30 I = 1,nmis
7530    ykmis(I) = SUM(c(NIYK(I),1:nz,S(1))*yk(ny+1:ny+nz))
76     +         + SUM(H(NIYK(I),1:nx,S(2))*STATE(1:nx))
77     +         + SUM(G(NIYK(I),1:nu,S(3))*U(1:nu))
78
79      DEALLOCATE (R,c,H,G,a,F)
80
81	RETURN
82	END
83