1 | c.......................................................................
   2 | c
   3 | c programme de traitement des donnees psi --> densites et bilans
   4 | c
   5 | c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
   6 | c Modele approche de rayonnement: methode pertubationnelle
   7 | c dst 1er ordre en kgaz et ksuies
   8 | c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
   9 | c
  10 | c
  11 | c.......................................................................
  12 | 
  13 |       program mosim
  14 | 
  15 | 
  16 |       implicit none
  17 | c     utilisation des variables du common, pas de passage de var
  18 |       include 'cecile.inc'
  19 | c pour def n_gaz_mx il faut propradia.inc
  20 |       include 'propradia.inc'
  21 |       include 'entre.inc'
  22 |       integer itmp1, i,j,m
  23 |       double precision bil ,var_bil, psid, supp
  24 |       double precision fvnew, fmnew, tempnew
  25 | 
  26 |       dimension fvnew(1:n_mx)
  27 |       dimension tempnew(0:n_mx+1)
  28 |       dimension fmnew(1:ngaz_mx, 0:n_mx+1)
  29 | 
  30 |       dimension bil(0:n_mx+1), var_bil(0:n_mx+1), psid(0:n_mx,0:n_mx)
  31 |       dimension supp(0:n_mx+1,0:n_mx+1)
  32 | 
  33 |       character*20 chemin
  34 | 
  35 |       print *, 'CHANGER LE CHEMIN DE LA SOL DE REFERENCE'
  36 |       chemin = './resultats/testref'
  37 | c la description des fichiers de sortie est dans cecile.out
  38 | c      print *, 'SORTIES fichiers: cecile.out'
  39 | c      print *
  40 | c      open (50,file='cecile.out',status='unknown')
  41 | c      write(50,*) 'Fichier de description des fichiers sorties'
  42 | c      write(50,*)
  43 | c      write(50,*) 'cecile.ecr: redirection de l''affichage ecran si b.'
  44 | c      write(50,*)
  45 | c      write(50,*) 'sorties ascii --> .out '
  46 | c      write(50,*) 'ind.out: Indices imax=n jmax=n kmax=n'
  47 | c      write(50,*) 'psi.out: Puissance nette echangees en W --- 
  48 | c     & Ecart type pne'
  49 | c      write(50,*) 'dpsidk.out: DERIV/kgaz --- ecart --- DERIV /ksuie --- 
  50 | c     & ecart'
  51 | c      write(50,*) 
  52 | c      write(50,*) 'sortie binaires  --> .oub'
  53 | c      close(50)
  54 | 
  55 | 
  56 | 
  57 | c ***********************************************************************
  58 | c ***********************************************************************
  59 | c Sorties ascii: ACCES SEQUENTIEL
  60 | c Leur format (commentaires ...)  doit s'adapter a celui utilise par le 
  61 | c logiciel de relecture.
  62 | c Faire une boucle sur chaque indice sauf si tableau mieux
  63 | 
  64 | 
  65 | 
  66 | 
  67 | c*************************************************************************
  68 | c mise a zero des variables
  69 |       
  70 |       do i=0,n+1
  71 |          do j=0,n+1
  72 | c            do r=1,n
  73 |                supp(i,j)=0.D+0
  74 | c            enddo
  75 |          enddo
  76 |       enddo
  77 | 
  78 | 
  79 | 
  80 | c ************************************************************************
  81 | c LECTURE FORMAT GNUPLOT <=> sur une ligne: param,..,function (un param varie/bloc)
  82 |       open(20,file=chemin // '/ind.out', status
  83 |      & ='unknown')
  84 |       read(20,*) n, rsup
  85 |       read(20,*) n, rsup
  86 |       read(20,*) n, rsup
  87 |       close(20)
  88 |  
  89 | c..................................................................................
  90 | c     open(20,file='psi.out',status='unknown', recl=21*(6))
  91 |       open(20,file=chemin // '/psi.out',status
  92 |      & ='unknown')
  93 | c     write(20,*) 'psi --- ecart type'
  94 |       do i = 0,n+1
  95 |               do j=0,n+1
  96 | c                 read(20,'(i4,i4,4(e22.13))') in, iin, r(in), r(iin)
  97 |                  read(20,*) in, iin, r(i), r(j)
  98 |      & ,psi(i,j), var_psi(i,j)
  99 |               enddo
 100 |               read(20,*)
 101 |       enddo
 102 |       close (20)
 103 | 
 104 | c...................................................................................
 105 | c     open(20,file='dpsidk.out',status='unknown',recl=21*(7))
 106 |       open(20,file=chemin // '/dpsi.out',status
 107 |      & ='unknown')
 108 | c     write(20,*)'DERIv psi/kgaz --- Ecart type --- DERIV psi/ksuie ---'
 109 |       do i=0,n+1
 110 | 		do j=0,n+1
 111 | 			do m=1,n
 112 | 				read(20,'(3(i4),4(e22.13))') in, iin, 
 113 |      &  itmp1, dpsi6dfv(i,j,m)
 114 | cccc     & , var_dpsi6dfv(i,j,m)
 115 |                                 !print *, dpsi6dfv(i,j,m)
 116 | 			enddo
 117 |                         read(20,*)
 118 | 		enddo
 119 |       enddo
 120 |       close(20)
 121 | c.....................................................................................
 122 | c     le fichier dpsidfm a ete rajoute le 14 novembre 1999
 123 |        open(20,file=chemin // '/dpsidfm.out',
 124 |      & status='unknown')
 125 |         do in=0,n+1
 126 | 		do iin=0,n+1
 127 | 			do itmp1=1,n
 128 | c				write(20,'(3(i4),2(e22.13))') in, iin,
 129 |                                 read(20,'(6(e22.13))') 
 130 |      &   dpsi6dfmco(in,iin,itmp1), 
 131 | cccc     &   var_dpsi6dfmco(in,iin,itmp1),
 132 |      &   dpsi6dfmco2(in,iin,itmp1), 
 133 | cccc     &   var_dpsi6dfmco2(in,iin,itmp1),
 134 |      &   dpsi6dfmh2o(in,iin,itmp1)
 135 | cccc     & , var_dpsi6dfmh2o(in,iin,itmp1)
 136 | 			enddo
 137 | 		enddo
 138 |       enddo
 139 |        close(20)
 140 | 
 141 | c.....................................................................................
 142 | c      open(20,file='dpsi.out',status='unknown',recl=21*(7))
 143 |       open(20,file=chemin // '/dpsidt.out',status
 144 |      & ='unknown')
 145 | c      write(20,*)'DERIv psi/kgaz --- Ecart type --- DERIV psi/ksuie ---'
 146 |       do i=0,n+1
 147 | 		do j=0,n+1
 148 | 			do m=1,n
 149 | c				read(20,'(3(i4),2(e22.13))') in, iin,
 150 |                                 read(20,*) in, iin,
 151 | c                           	write(20,*) in, iin,
 152 |      &   itmp1, dpsi6dt(i,j,m)
 153 | cccc     & , var_dpsi6dt(i,j,m)
 154 | 			enddo
 155 | c simple retour a la ligne 
 156 |                         read(20,*)
 157 | 		enddo
 158 |       enddo
 159 |       close(20)
 160 | 
 161 | 
 162 | 
 163 | 
 164 | c-----------------------------------------------------------------
 165 | c lectures des variables dans l'etat de refenrence
 166 | C LECTURE DE LA DIFFERENCE dfv
 167 |        open(20,file=chemin // '/fv.in',status='unknown')
 168 |        do i=1,n 
 169 |           read (20,*) fv(i)
 170 |        enddo
 171 | c
 172 |        close(20)
 173 | c les fractions molaires
 174 |         open(20,file=chemin // '/fmco.in',
 175 |      &     status='unknown')
 176 |        do i=1,n 
 177 |           read (20,*) fm(1,i)
 178 |        enddo
 179 | c
 180 |        close(20)
 181 |         open(20,file=chemin // '/fmco2.in',
 182 |      &      status='unknown')
 183 |        do i=1,n 
 184 |           read (20,*) fm(2,i)
 185 |        enddo
 186 | c
 187 |        close(20)
 188 |         open(20,file=chemin // '/fmh2o.in',
 189 |      &      status='unknown')
 190 |        do i=1,n 
 191 |           read (20,*) fm(3,i)
 192 |           !print *, fm(3,i)
 193 |        enddo
 194 | c
 195 |        close(20)
 196 | C LECTURE DE LA temp
 197 |        open(20,file=chemin // '/temperat.in',
 198 |      &     status='unknown')
 199 |        do i=0,n+1 
 200 |           read (20,*) temp(i)
 201 |        enddo
 202 | 
 203 |        close(20)
 204 | c---------------------------------------------------------------------
 205 | c lecture nouvaeu profil qui un jour sera passer par egolfopoulos
 206 | C LECTURE DE LA DIFFERENCE dfv
 207 |        open(20,file='./fv.in',status='unknown')
 208 |        do i=1,n 
 209 |           read (20,*) fvnew(i)
 210 |        enddo
 211 | c
 212 |        close(20)
 213 | c les fractions molaires
 214 |         open(20,file='./fmco.in',status='unknown')
 215 |        do i=1,n 
 216 |           read (20,*) fmnew(1,i)
 217 |        enddo
 218 | c
 219 |        close(20)
 220 |         open(20,file='./fmco2.in',status='unknown')
 221 |        do i=1,n 
 222 |           read (20,*) fmnew(2,i)
 223 |        enddo
 224 | c
 225 |        close(20)
 226 |         open(20,file='./fmh2o.in',status='unknown')
 227 |        do i=1,n 
 228 |           read (20,*) fmnew(3,i)
 229 |           !print *, fmnew(3,i), (fmnew(3,i)-fm(3,i))
 230 |        enddo
 231 | c
 232 |        close(20)
 233 | C LECTURE DE LA temp
 234 |        open(20,file='./temperat.in',status='unknown')
 235 |        do i=0,n+1 
 236 |           read (20,*) tempnew(i)
 237 |        enddo
 238 | 
 239 |        close(20)
 240 | 
 241 |     
 242 | c......................................................................
 243 | c sortie binaire: ACCES DIRECT (plus universel)
 244 | c cas ou les sorties ascii ne sont pas bien mise en forme et ou il faut
 245 | c tout lire rapidement
 246 | c chaque ligne est un enregistrement qu'on peut lire d'un bloc
 247 | c enregistrement depend de la machine: sun double precision ici
 248 | c on peut lire les lignes dans le desordre
 249 |      
 250 | c numero ligne ou on enregistre moins un 
 251 | ccc      itmp5= 0
 252 | ccc      open(20,file='psi.oub',access='direct',recl=8 
 253 | ccc     &*(n+1+1))
 254 | ccc     write(20,*) 'Matrice des psi'
 255 | ccc      do in=0,n+1
 256 | ccc	      write(20,rec=in+1)  (psi(in,iin), iin=0,n+1)
 257 | ccc      enddo 
 258 | c derniere valeur de rec
 259 | ccc      itmp5=in+1
 260 | 
 261 | c     write(20,*)'Matrice des ecarttypes de psi')
 262 | ccc      do in=0+itmp5,n+1+itmp5
 263 | ccc	write(20,rec=in+1) (var_psi(in,iin), iin=0,n+1)
 264 | ccc      enddo
 265 | 
 266 | c
 267 | c gerer bloc acces direct connitre numero de ligne
 268 | c
 269 | ccc      close(20)
 270 | 
 271 | 
 272 | 
 273 | c ************************************************************************
 274 | c TRAITEMENT PAR LE DL
 275 | 
 276 | 
 277 | 
 278 |        do i=0,n+1 
 279 |           do j=0,n+1
 280 | 
 281 | c         print *,'Afficaghage',psi(i,i), dpsi6dfv(i,i,i), dfv(i)
 282 |              do m=1,n
 283 | c                supp(i,j)=supp(i,j)+dpsi6dfv(i,j,m)*(fvnew(m)- fv(m))
 284 |                 supp(i,j)=supp(i,j)+dpsi6dfv(i,j,m)*(fvnew(m)- fv(m))
 285 |      &                         +dpsi6dfmco(i,j,m)*(fmnew(1,m)-fm(1,m))
 286 |      &                         +dpsi6dfmco2(i,j,m)*(fmnew(2,m)-fm(2,m))
 287 |      &                         +dpsi6dfmh2o(i,j,m)*(fmnew(3,m)-fm(3,m))
 288 |      &                         +dpsi6dt(i,j,m)*(tempnew(m)-temp(m))
 289 |              enddo
 290 |              !if( i.eq.3) then
 291 |                 
 292 |              !print *, 'i, j, sup:  ', i, j, supp(i,j)
 293 |              !read *
 294 |              !endif
 295 |              psi(i,j)=psi(i,j)+supp(i,j)
 296 |              
 297 |           enddo
 298 | 
 299 | cc          print *, (fvnew(i)-fv(i))
 300 | cc          print *, (fmnew(3,i)-fm(3,i))
 301 | cc          print *, (fmnew(2,i)-fm(2,i))
 302 |           print *, (fmnew(1,i)-fm(1,i))
 303 | cc          print *, (tempnew(i)-temp(i))
 304 |       enddo
 305 | 
 306 | 
 307 | 
 308 | 
 309 |       call fich
 310 | 
 311 | 
 312 |       end
 313 | 
 314 |