source: LMDZ.3.3/trunk/libf/phylmd/testphys1d.F @ 346

Last change on this file since 346 was 248, checked in by lmdz, 23 years ago

Routines pour le 1D F.Hourdin
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 18.4 KB
Line 
1      PROGRAM gcm1d
2      IMPLICIT NONE
3
4#include "dimensions.h"
5#include "dimphy.h"
6#include "YOMCST.h"
7#include "comg1d.h"
8#include "clesphys.h"
9#include "control.h"
10
11c    Arguments :
12c    -----------
13
14      integer ngrid,nlayer,longcles,nqmax
15      parameter (longcles=20,nqmax=3)
16cfleur on  ne se sert pas de nqmax mais de nqmx defini dnas dimensions.h depend de la
17c      compilation
18      integer radpas
19      integer nday,it,iflag_con,unit,i,itap,n_cooling
20      integer iphys_ver,iadv_tvl,i_cvg,i_hum
21
22      real fnday
23      real h,kappa,z1,z2,sbid,sigbid
24      real plev(klev+1),play(klev), psol
25      real temp(klev),u(klev),v(klev),tsurf,q(klev,nqmx),w(klev+1)
26      real ug(klev),vg(klev)
27      real du(klev),dv(klev),dt(klev),dpsrf,dq(klev,nqmx)
28      real du_dyn(klev),dv_dyn(klev)
29     :    ,dt_dyn(klev),dq_dyn(klev,nqmax)
30      real phis,presnivs(klev),clesphy0(longcles)
31      real time,timestep,ecritphy,day,tho
32      real co2_ppm,solaire
33      real rlat,rlon,tsol,radsol,psol_f,tsol_f,qsol_f,qsol
34      real rugmer,rugsrel,snow,agesno,deltat,zmea,zstd,zsig,sn
35      real zgam,zthe,zpic,zval
36      real ema_sig(klev),ema_w(klev),ncst_cbmf,ema_cbmfz,ema_pcb
37     : ,zz_f(klev),vu_f(klev),vv_f(klev),t_f(klev),q_f(klev,nqmax)
38
39      real temp0(klev),q0(klev,3)
40
41      real dt_cooling(klev),dq_cooling(klev)
42      real d_t_cool(klev),d_q_cool(klev)
43      real d_t_adv(klev),d_q_adv(klev,nqmax)
44      real d_t_cvg(klev),d_q_cvg(klev)
45      real ht(100),hq(100),hw(100)
46
47      real phy_nat(360)
48      real phy_alb(360)
49      real phy_sst(360)
50      real phy_bil(360)
51      real phy_rug(360)
52      real phy_ice(360)
53
54      logical cycle_diurne,soil_model,new_oliq,ok_orodr
55     :       ,ok_orolf,ok_limitvrai
56
57      logical firstcall,lastcall,flag_cool
58      logical itsourcecont
59c      itsourcecont permet de choisir entre les sources pour plus de 5 traceurs
60      character*80 ans,file_forctl, file_fordat, file_start,file
61      character*2 cnbl
62
63c-----------------------------------------------------------------------
64
65      COMMON/comvert/
66     * s(llm),sig(llm+1),ds(llm),dsig(llm),dsig1(llm),sdsig(llm)
67     . ,sig_s(llm)
68
69      REAL s,sig,ds,dsig,dsig1,sdsig,sig_s
70
71c-----------------------------------------------------------------------
72c INCLUDE 'temps.h':
73
74      COMMON/temps/itaufin,dtd,
75     s  day_ini,day_end,anne_ini
76
77      INTEGER  itaufin
78      INTEGER*4 day_ini,day_end,anne_ini
79      REAL     dtd
80
81c-----------------------------------------------------------------------
82c   dynamical tendencies
83
84      INTEGER l,ierr,aslun,nlevel,iq,ll
85
86      REAL longitude,latitude
87      REAL zlay(klev),phi(klev)
88      REAL paire
89
90      DATA latitude,longitude/0.,0./
91
92c-----------------------------------------------------------------------
93c    Initialisations  des constantes
94c    -------------------------------
95
96c constantes
97c ----------
98
99      time=0.
100      tho=3600.
101      it=0
102
103      call suphec
104
105c parametres lus dans execution_1D:
106c ---------------------------------
107
108      read(*,*) fnday
109      print*, 'fnday',fnday
110      read(*,*) ecritphy
111      print*, 'ecritphy',ecritphy
112      read(*,*) timestep
113      print*, 'timestep',timestep
114      read(*,*) iflag_con
115      print*, 'iflag_con',iflag_con
116      read(*,*) flag_cool
117      print*, 'flag_cool',flag_cool
118      read(*,'(a)') cnbl
119      print*, 'cnbl',cnbl
120      read(*,*) radpas   
121      print*, 'radpas',radpas
122      read(*,*) anne_ini ! sans importance mais il faut qqchose
123      print*, 'anne_ini',anne_ini
124     
125c si ecritphy=0: on ecrit la physique a chaque pas de temps:
126
127      if (ecritphy.eq.0.) then
128       ecritphy = timestep/rday
129      endif
130      write(*,*) 'ECRITPHY = ',ecritphy
131
132      firstcall=.true.
133      lastcall=.false.
134      day=1.
135      itsourcecont=.true.
136
137      open(78,file='transil.dat',form='unformatted',
138     s    access='direct',recl=4*llm*llm)
139
140      ngrid=1
141      IF (ngrid.NE.klon) THEN
142         PRINT*,'STOP in inifis'
143         PRINT*,'Probleme de dimenesions :'
144         PRINT*,'ngrid     = ',ngrid
145         PRINT*,'klon   = ',klon
146         STOP
147      ENDIF
148      nlayer=klev
149      nlevel=nlayer+1
150
151c   -------------------------------------------------------------------
152c   Initialisation de la discretisation verticale
153c   ---------------------------------------------
154c sb      CALL disvert(llm,rkappa,sig,dsig,s,ds,dsig1,sdsig)
155
156c   -------------------------------------------------------------------
157c   Profils initiaux.
158c   -----------------
159
160      unit = 80
161      open(unit,file='start'//cnbl//'.data',form='formatted')
162      read(unit,*)
163      read(unit,*)
164      read(unit,*) (play(l),l=1,nlayer)
165      read(unit,*) (plev(l),l=1,nlevel)
166      read(unit,*) (temp(l),l=1,nlayer)
167      read(unit,*) (q(l,1),l=1,nlayer)
168      read(unit,*) (q(l,2),l=1,nlayer)
169      read(unit,*) (ema_sig(l),l=1,nlayer)
170      read(unit,*) (ema_w(l),l=1,nlayer) ! reinitialise + loin???
171      read(unit,*) ncst_cbmf, ema_cbmfz, ema_pcb
172      read(unit,*) (zz_f(l),l=1,nlayer)
173      read(unit,*) (vu_f(l),l=1,nlayer)
174      read(unit,*) (vv_f(l),l=1,nlayer)
175      close(unit)
176      write(*,*) 'Lecture fichier start.data ok'
177
178      DO l = 1, nlayer
179         u(l) = vu_f(l)
180         v(l) = vv_f(l)
181c##################
182cfleur ATTENTION q(3) ce nest plus la glace mais le radon si on active phytrac
183c##################3
184         q(l,3) = 0. ! on initialise la glace a zero
185         PRINT*,plev(l),play(l)
186c on recopie le profil initial dans les champs de forcing:
187         t_f(l) = temp(l)
188         q_f(l,1) = q(l,1)
189         q_f(l,2) = q(l,2)
190         q_f(l,3) = q(l,3)
191      ENDDO
192
193
194c Lecture/creation des conditions aux limites:
195c --------------------------------------------
196
197      unit = 81
198      open(unit,file='startphy.data',form='formatted')
199      read(unit,*)
200      read(unit,*) co2_ppm
201      read(unit,*) solaire
202      read(unit,*) rlat
203      read(unit,*) tsol
204      read(unit,*) radsol
205      close(unit)
206      write(*,*) 'Lecture fichier startphy.data ok'
207
208      unit = 82
209      open(unit,file='condsurf.data',form='formatted')
210      read(unit,*) psol_f
211      read(unit,*) tsol_f
212      read(unit,*) qsol_f
213      close(unit)
214      write(*,*) 'Lecture fichier condsurf.data ok'
215
216      if (play(1).lt.10000.) then
217        do l = 1, nlayer
218         play(l) = play(l)*100.
219        enddo
220      endif
221      if (plev(1).lt.10000.) then
222        do l = 1, nlevel
223         plev(l) = plev(l)*100.
224        enddo
225      endif
226
227      if (abs(psol_f-plev(1)) .gt. 1.) then
228         print *,' Incompatibilite entre psol et profil'
229     :   ,' de pression'
230         print *,' psol = ',psol_f,' plev(1) = ',plev(1)
231         stop
232      endif
233
234      tsol    = tsol_f ! temp au sol prise dans condsurf.data
235      psol    = psol_f ! pression au sol prise dans condsurf.data
236      qsol    = qsol_f
237      rugmer  = 0.0001 ! valeur de cchlim.data
238      rugsrel = 0.0    ! (rugsrel = rugoro)
239      snow    = 0.0
240      sn=0.
241      agesno  = 50.0
242      rlon    = 0.0
243      deltat  = 0.0    ! ne sert que pour les slab_ocean
244      phis = 0.
245      zmea = 0.
246      zstd = 0.
247      zsig = 0.
248      zgam = 0.
249      zthe = 0.
250      zpic = 0.
251      zval = 0.
252
253      do i=1,360
254      phy_sst(i) = tsol
255      phy_nat(i) = 0.0  ! 0=ocean libre, 1=land, 2=glacier, 3=banquise
256      phy_alb(i) = 0.15 ! albedo land only (old value condsurf_jyg=0.3)
257      phy_bil(i) = 1.0  ! ne sert que pour les slab_ocean
258      phy_rug(i) = 0.1  ! longueur rugosite utilisee sur land only
259      phy_ice(i) = 0.0  ! fraction de glace (?)
260      enddo
261
262      call writelim
263     s   (phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice)
264
265
266
267c controles du run (en 3D: lus dans run.def) :
268
269      cycle_diurne = .FALSE.
270      soil_model   = .FALSE.
271      new_oliq     = .FALSE.
272      ok_orodr     = .FALSE.
273      ok_orolf     = .FALSE.
274      ok_limitvrai = .FALSE.
275
276      do i = 1, longcles
277       clesphy0(i) = 0.
278      enddo
279
280      nbapp_rad = NINT(86400./radpas/timestep)
281
282      clesphy0(1) = FLOAT( iflag_con )
283      clesphy0(2) = FLOAT( nbapp_rad )
284      if (cycle_diurne ) clesphy0(3) = 1.
285      if (  soil_model ) clesphy0(4) = 1.
286      if (    new_oliq ) clesphy0(5) = 1.
287      if (    ok_orodr ) clesphy0(6) = 1.
288      if (    ok_orolf ) clesphy0(7) = 1.
289      if ( ok_limitvrai) clesphy0(8) = 1.
290
291
292c -------------------------------------------------------------------
293c Discretisation verticale:
294c -------------------------
295
296c calcul sig,dsig,s,ds,dsig1,sdsig a partir des play, plev lus:
297
298      do l = 1, nlevel
299       sig(l)= plev(l)/plev(1)
300      enddo
301
302      do l = 1, nlayer
303       s(l)= sig(l)**RKAPPA
304      enddo
305
306      do l = 2, nlayer
307       ds(l) = s(l-1) - s(l)
308      enddo
309      ds(1)  = 1. - s(1)
310
311      do l = 1, nlayer
312       dsig(l)  = sig(l)-sig(l+1)
313       sdsig(l) = s(l) * dsig(l)
314       dsig1(l) = 1./dsig(l)
315       presnivs(l) = play(l)/100./100.
316      enddo
317
318      print*,'Diagnostique de la discretisation verticale'
319
320      print*
321      h=7.
322      kappa = RKAPPA
323      print*,'comparaison de sig(l) et (s(l)+s(l+1))/2)**(1/K)'
324      do 14 l=1,llm-1
325         sigbid=(0.5*(s(l)+s(l+1)))**(1./kappa)
326         print*,'sig(',l+1,')  = ',sig(l+1),
327     :           '    valeur approchee :',sigbid,'   ',dsig(l)
32814    continue
329      print*
330      print*,'comparaison de s(l) et (sig(l)+sig(l+1))/2)**K'
331      do 15 l=1,llm
332         sbid=(0.5*(sig(l+1)+sig(l)))**kappa
333         print*,'  s(',l,')  = ',s(l),
334     :           '    valeur approchee :',sbid
33515    continue
336
337      print*,'Altitude approchee z,dz'
338      print*
339      z1=0.
340      print*,'   l       Z      DZ      Ztop   dsig'
341      DO 18 l=1,llm-1
342         z2=-h*log(sig(l+1))
343         write(*,'(i5,3x,4f8.4)') l,-h*log(s(l))/kappa,z2-z1,z2
344     :    ,dsig(l)
345         write(14,'(3x,i5,1f10.4)') l,-h*log(s(l))/kappa
346         z1=z2
34718    CONTINUE
348      write(*,'(i5,3x,3f8.4)') l,-h*log(s(llm))/kappa
349      write(14,'(3x,i5,1f10.4)') l,-h*log(s(llm))/kappa
350
351c      print*,'Correspondance approx pression-altitude:'
352c      z1 = 0.
353c      do l = 1, llm
354c       z1 = z1 + 287.*temp(l)/9.81/play(l)*(plev(l)-plev(l+1))
355c       write(*,*) l,play(l),z1
356c      enddo
357
358c -------------------------------------------------------------------
359c Ecriture de l'etat initial:
360c ---------------------------
361
362c     call ini_fis(timestep,radpas,iflag_con,anne_ini
363c    :            ,co2_ppm,solaire,rlat,rlon,tsol,deltat
364c    :            ,qsol,snow,radsol,rugmer,agesno,zmea,zstd,zsig,zgam
365c    :            ,zthe,zpic,zval,rugsrel
366c    :            ,phy_sst,phy_nat,phy_alb,phy_bil,phy_rug,phy_ice)
367
368      call physdem(rlon, rlat, timestep,radpas,co2_ppm,
369     .                   solaire,tsol, qsol,
370     .                   sn, radsol, deltat, rugmer,
371     .                   agesno, zmea, zstd, zsig,
372     .                   zgam, zthe, zpic, zval,
373     .                   rugsrel)
374
375
376
377c -------------------------------------------------------------------
378c Options de la simulation:
379c -------------------------
380
381      unit = 83
382      open(unit,file='version.data',form='formatted')
383      read(unit,*) iphys_ver,iadv_tvl,i_cvg,i_hum
384      close(unit)
385      write(*,*) 'Lecture fichier version.data ok'
386
387      write(*,*) 'physiq # ',iphys_ver,' (0=std; 1=forced glob cis;
388     : 2=forced loc cis)'
389      write(*,*) 'advection switch : ',iadv_tvl,' (0=no adv ;
390     : 1=forced adv)'
391      write(*,*) 'convergence switch : ',i_cvg,' (0=no conv;
392     : 1=forced conv)'
393
394
395
396c Eventuellement, preparation de la "bulle froide":
397c -------------------------------------------------
398
399      IF (flag_cool) THEN
400      unit = 84
401      open(unit,file='cool_buble.data',form='formatted')
402      read(unit,*) n_cooling
403      do l = 1, nlayer
404      read(unit,*) dt_cooling(l),dq_cooling(l)
405      enddo
406      close(unit)
407      write(*,*) 'Lecture fichier cool_buble.data ok'
408      ENDIF
409
410
411c Eventuellement, preparation du forcage par la convergence:
412c ----------------------------------------------------------
413
414      IF (i_cvg .EQ. 1) THEN
415      file_forctl = 'forcing.ctl'
416      file_fordat = 'forcing.dat'
417      file_start  = 'start'//cnbl//'.data'
418
419      call copie(klev,play,psol,file_forctl)
420      call get_uvd2(itap,file_forctl,file_fordat,file_start
421     :             ,ht,hq,hw)
422      ENDIF
423
424
425c-----------------------------------------------------------------------
426c  initialisation pour GRADS-1D
427c  ----------------------------
428
429      g1d_nlayer=nlayer
430      g1d_nomfich='grads1d.dat'
431      g1d_unitfich=30
432      g1d_nomctl='grads1d.ctl'
433      g1d_unitctl=31
434      g1d_premier=.true.
435
436      file='sort'
437      call inigrads(1,1
438     s  ,0.,1.,-2.,2.,1,0.,-2.,2.,1.
439     s  ,llm,presnivs,1000.
440     s  ,1800.,file,'Diagconvect')
441      print*,'Fin de Initialisation de wrgras'
442
443
444c=======================================================================
445c   DEBUT DE L'INTEGRATION TEMPORELLE:
446c   ==================================
447
448      itap = 1
449
4501     continue
451
452c  calcul du geopotentiel:
453c  -----------------------
454
455      phi(1)=RD*temp(1)*
456     :(plev(1)-play(1))/(.5*(plev(1)+play(1)))
457      do l = 1, nlayer-1
458         phi(l+1)=phi(l)+RD*(temp(l)+temp(l+1))*
459     :   (play(l)-play(l+1))/(play(l)+play(l+1))
460      enddo
461
462c        PRINT*,'altitude km  et T (K)'
463c        WRITE(6,'(2f10.2)') (.001*phi(l)/g,temp(l),l=1,nlayer)
464
465
466c pour etre coherent avec LMDZ.3, on passe 2 arguments
467c en plus a physiq: phis et aire (pas utilises en 1D):
468      paire = 1. ! aire de la maille
469
470
471c  eventuellement, induction de la convection par chauffage du
472c  bas et refroidissement du haut de la couche limite:
473c  ------------------------------------------------------------
474
475      IF (flag_cool) THEN
476
477      call cool_pool(itap
478     e           ,n_cooling,dt_cooling,dq_cooling
479     s           ,d_t_cool,d_q_cool)
480
481      do l = 1, nlayer
482       temp(l) = temp(l) + d_t_cool(l)
483       q(l,1)  = q(l,1)  + d_q_cool(l)
484      enddo
485
486      ELSE
487
488      do l = 1, nlayer
489       d_t_cool(l) = 0.
490       d_q_cool(l) = 0.
491      enddo
492
493      ENDIF
494
495
496c  eventuellement, ajouter une tendance relative a une
497c  translation du domaine (ex: pour suivre une ligne de grains):
498c  -------------------------------------------------------------
499
500      IF (iadv_tvl .EQ. 1) THEN
501
502      call advect_tvl(timestep,temp,q,vu_f,vv_f,t_f,q_f
503     :               ,d_t_adv,d_q_adv)
504
505      do l = 1, nlayer
506       temp(l) = temp(l) + d_t_adv(l)
507       q(l,1)  = q(l,1)  + d_q_adv(l,1)
508       q(l,2)  = q(l,2)  + d_q_adv(l,2)
509       q(l,3)  = q(l,3)  + d_q_adv(l,3)
510      enddo
511
512      ELSE
513
514      do l = 1, nlayer
515       d_t_adv(l) = 0.
516       d_q_adv(l,1) = 0.
517       d_q_adv(l,2) = 0.
518       d_q_adv(l,3) = 0.
519      enddo
520
521      ENDIF
522
523
524c  eventuellement, ajouter une tendance relative a la
525c  convergence de grande echelle:
526c  -------------------------------------------------------------
527
528      IF (i_cvg .EQ. 1) THEN
529
530      call get_uvd(itap,timestep,tsol,qsol,file_fordat,ht,hq,hw)
531
532      do l = 1, nlayer
533
534       temp0(l) = temp(l) ! memoire de "l'avant dynamique"
535       q0(l,1) = q(l,1)
536       q0(l,2) = q(l,2)
537       q0(l,3) = q(l,3)
538
539       d_t_cvg(l) = ht(l) * timestep
540       d_q_cvg(l) = hq(l) * timestep
541       temp(l) = temp(l) + d_t_cvg(l)
542       q(l,1)  = q(l,1)  + d_q_cvg(l)
543c       write(*,*)'d_t_cvg,d_q_cvg:',d_t_cvg(l),d_q_cvg(l)
544       if (q(l,1).lt.0.) then ! sb
545       print*,'OVAP negative dans main!'
546       print*,'itap,l,q,d_q_cvg:',itap,l,q(l,1),d_q_cvg(l),hq(l)
547       q(l,1)  = MAX(q(l,1),1.e-10) ! evite les humidites negatives
548       endif
549      enddo
550
551      if (itap.ge.30) then
552        print*,'hq(l),ht(l),dq,dt,q:'
553     :  ,hq(16),ht(16),d_q_cvg(16),d_t_cvg(16),q(16,1)
554      endif
555
556      ELSE
557
558      do l = 1, nlayer
559       d_t_cvg(l) = 0.
560       d_q_cvg(l) = 0.
561      enddo
562
563      ENDIF
564
565      do l = 1, nlayer
566       dt_dyn(l) = d_t_cvg(l) / timestep
567       dq_dyn(l,1) = d_q_cvg(l) / timestep
568       dq_dyn(l,2) = 0.0
569       dq_dyn(l,3) = 0.0
570      enddo
571
572
573c   calcul des tendances physiques:
574c   -------------------------------
575
576      if (itsourcecont) then
577        do iq=5,nqmx
578          ll=min(iq-4,18)
579          if(firstcall) then
580c          q(ll,iq)=5/((plev(ll)-plev(ll+1))/rg)
581c           q(ll,iq)=5
582c          q(ll,iq)=5*play(ll)/((plev(ll)-plev(ll+1))/rg)/(RD*temp(ll))
583            q(ll,iq)=5*play(ll)/(RD*temp(ll))
584          endif
585        print*,'q ll iq',q(ll,iq),ll,iq
586        enddo
587      else
588        do iq=5,nqmx
589          ll=min(iq-4,18)
590          q(ll,iq)=q(ll,iq)+1
591             do l=1,llm
592               q(l,iq)=q(l,iq)*(1-timestep/tho)
593        print*,'q ll iq',q(ll,iq),ll,iq
594             enddo   
595        enddo
596 
597      endif
598
599c      CALL physiq(ecritphy,ngrid,nlayer,nqmx,
600c     :            firstcall,lastcall,
601c     :            day,day,time,timestep,
602c     :            plev,play,phi,phis,paire,
603c     :            u,v,temp,q,
604c     :            du_dyn, dv_dyn,  dt_dyn, dq_dyn,
605c     :            w,
606c     :            du,dv,dt,dq,dpsrf)
607
608c     CALL physiq (ecritphy,ngrid,nlayer,nqmx,firstcall,lastcall,
609c    ,   day,day,time,timestep,plev,play,phi,phis,paire,
610c     ,             presnivs,clesphy0,u,v,temp,q, 
611c pour calculer proprement la tendance de cape liee a la
612c dynamique, il faut entrer aussi temp0 et q0:
613c    ,             presnivs,clesphy0,u,v,temp,q,temp0,q0, 
614c    ,             du_dyn,dv_dyn,dt_dyn,dq_dyn,w,
615c - sorties
616c    s             du,dv,dt,dq,dpsrf)
617
618      print*,'PAS DE TEMPS ',timestep
619      print*,'ATTENTION!!! Il faudra passer temp0 et q0'
620      CALL physiq(ngrid,nlayer,nqmx,
621     :            firstcall,lastcall,
622     :            day,day,time,timestep,
623     :            plev,play,phi,phis,paire,presnivs,clesphy0,
624     :            u,v,temp,q,
625     :            w,
626     :            du,dv,dt,dq,dpsrf)
627
628
629
630
631      firstcall=.false.
632
633
634c   Ajout des tendances
635c   -------------------
636      DO l=1,nlayer
637         u(l)=u(l)+timestep*du(l)
638         v(l)=v(l)+timestep*dv(l)
639
640         temp(l)=temp(l)+timestep*dt(l)
641      ENDDO
642
643      do iq=1,nqmx
644         do l=1,nlayer
645            q(l,iq)=q(l,iq)+timestep*dq(l,iq)
646         enddo
647      enddo
648c      write(78,rec=itap) ((q(l,iq),l=1,llm),iq=5,22)
649      write(78,rec=itap) ((q(l,iq),iq=5,22),l=1,llm)
650      itap = itap + 1
651
652
653         time=time+timestep/rday
654         it=it+1
655         if(mod(it,2000).eq.0) print*,'TIME=',time
656         if(time.gt.fnday) then
657            CALL endg1d(1,nlayer,play,int(time/ecritphy),timestep)
658            stop
659         endif
660
661         GOTO 1
662
663      END
664
665c=======================================================================
666c=======================================================================
667c    FIN DU PROGRAMMES
668c    CI-DESSOUS, QUELQUES SOUS-PROGRAMMES UTILS
669c=======================================================================
670c=======================================================================
671
672#include "1DUTILS.h"
673#include "1Dconv.h"
Note: See TracBrowser for help on using the repository browser.