source: trunk/hftp/iasi/aerosol/prog/gads.f @ 1

Last change on this file since 1 was 1, checked in by cbipsl, 18 years ago

Geisa inital import
File size: 61.3 KB
Line 
1      program gads
2     
3ccccc ------------------------------------------------------------------c
4c     create global distributions of microphysical and optical aerosol  c
5c     properties on the base of the GADS database.                      c
6c                                                                       c
7c     version 2.0                                                       c
8c     -----------                                                       c
9c     11.02.97                                                          c
10c                                                                       c
11c     version 2.1                                                       c
12c     -----------                                                       c
13c     19.11.97 - new file format in ../optdat/                          c
14c              - dimensions for phase function now 112.                 c
15c                                                                       c
16c     version 2.2                                                       c
17c     -----------                                                       c
18c     21.01.98 - files winter.dat, summer.dat updated                   c
19c              - output improved                                        c
20c              - optical depth corrected                                c
21c                                                                       c
22c     version 2.2a                                                      c
23c     -----------                                                       c
24c     06.03.98 - error with calculation of mass corrected.              c
25c                                                                       c
26c                                                                       c
27c     06.03.98                                                 M. Hess  c
28ccccc ------------------------------------------------------------------c
29
30      character*1 om
31
32      print*,'********************************************************'
33      print*,'*             Global Aerosol Data Set 2.2a             *'
34      print*,'*  --------------------------------------------------  *'
35      print*,'*                                                      *'
36      print*,'*  GADS is described in:                               *'
37      print*,'*  P. Koepke, M. Hess, I. Schult,                      *'
38      print*,'*                  and E.P. Shettle (1997):            *'
39      print*,'*  "Global Aerosol Data Set"                           *'
40      print*,'*  submitted to Theoret. and Appl. Climate             *'     
41      print*,'*                                                      *'
42      print*,'*  preprint available as Report No. 243,               *'
43      print*,'*  Max-Planck-Institut fuer Meteorologie, Hamburg      *'
44      print*,'*                                                      *'
45      print*,'*  last update: 06.03.98                     M. Hess   *'   
46      print*,'********************************************************'
47      print*,' '
48      print*,' '
49           
50      print*,'  '
51      write(*,114)
52  114 format(' do you want to get (m)icrophysical or (o)ptical data?')
53      read (*,'(a)') om
54     
55      if (om.eq.'m') then
56         call mic
57      else
58         call opt
59      end if
60     
61
62      stop 'Normal end of GADS'
63      end
64     
65      subroutine mic
66
67ccccc -----------------------------------------------------------------c
68c     veraenderte Version von GLOPLO.FOR                               c
69c                                                                      c
70c     Es werden number densities und mix. rat. aus Rohdaten geplottet  c
71c                                                                      c
72c     Es sind maximal 5 Komponenten pro Typ vorgesehen.                c
73c                                                                      c
74c     Stand: 03.03.97                                         M. Hess  c
75ccccc -----------------------------------------------------------------c
76
77      real      w(73,37)
78      real      vn(10),dens(10),w2(73,37)
79      integer   x(73),y(37),yi,xi,rh
80      integer   t1,t2,t3,t4,t5
81
82      character*1  ws
83      character*2  at
84      character*5  parout(23)
85      character*7  season
86      character*9  file
87      character*31 parnam,parn(24)
88
89      logical ende
90
91      data parout /'tonde','inson','wason','sootn',
92     *             'seaan','seacn','minnn','minan','mincn',
93     *             'mintn','sulfn','aetyp','insom','wasom','sootm',
94     *             'seaam','seacm','minnm',
95     *             'minam','mincm','mintm','sulfm','prtyp'/
96
97      data parn /'number density (part./cm**3)  ',
98     *           'insoluble, (part./cm**3)      ',
99     *           'water soluble, (part./cm**3)  ',
100     *           'soot, (part./cm**3)           ',
101     *           'sea salt (acc.), (part./cm**3)',
102     *           'sea salt (coa.), (part./cm**3)',
103     *           'mineral (nuc.), (part./cm**3) ',
104     *           'mineral (acc.), (part./cm**3) ',
105     *           'mineral (coa.), (part./cm**3) ',
106     *           'mineral (tra.), (part./cm**3) ',
107     *           'sulfate, (part./cm**3)        ',
108     *           'Aerosol Types                  ',
109     *           'insoluble, mikrogr/(m**3)      ',
110     *           'water soluble, mikrogr/(m**3)  ',
111     *           'soot, mikrogr/(m**3)           ',
112     *           'sea salt (acc.), microgr/(m**3)',
113     *           'sea salt (coa.), microgr/(m**3)',
114     *           'mineral (nuc.), microgr/(m**3) ',
115     *           'mineral (acc.), microgr/(m**3) ',
116     *           'mineral (coa.), microgr/(m**3) ',
117     *           'mineral (tra.), microgr/(m**3) ',
118     *           'sulfate, microgr/(m**3)        ',
119     *           'Profile Typ                    ',
120     *           '                               '/
121
122      data vn   /1.19e7,7.43e2,5.98e1,3.64e5,
123     *           1.01e8,1.07e4,2.13e6,1.23e8,7.37e6,
124     *           1.34e4/
125      data dens /2.0,1.8,1.0,2.2,2.2,2.6,2.6,2.6,2.6,1.7/
126
127ccccc -----------------------------------------------------------------c
128c     Abfrage, was geplottet werden soll                               c
129ccccc -----------------------------------------------------------------c
130
131 1001 print*,'  '
132      write(*,154)
133  154 format(' (w)inter or (s)ummer? ')
134      read (*,'(a)') ws
135      if (ws.eq.'w') then
136         open(7,file='../glodat/winter.dat')
137         season='winter '
138      else if (ws.eq.'s') then
139         open(7,file='../glodat/summer.dat')
140         season='summer '
141      else
142         print*,' wrong input! try again!'
143         goto 1001
144      end if
145
146ccccc -----------------------------------------------------------------c
147c     Endlosschleife ueber restliches Programm                         c
148ccccc -----------------------------------------------------------------c
149
150      ende=.false.
151 1100 do while (.not.ende)
152
153      print*,'  '
154      print*,'  '
155      print*,' choose parameter to extract from database:'
156      print*,' '
157      print*,'               NUMBER DENSITIES:          MASS (rh=0%):'
158      print*,'  '
159      print*,'           (1) total number density'
160      print*,'    '
161      print*,'           (2) insoluble             (13) insoluble     '
162      print*,'           (3) watersoluble          (14) watersoluble  '
163      print*,'           (4) soot                  (15) soot          '
164      print*,'           (5) sea salt (acc)        (16) sea salt (acc)'
165      print*,'           (6) sea salt (coa)        (17) sea salt (coa)'
166      print*,'           (7) mineral (nuc)         (18) mineral (nuc) '
167      print*,'           (8) mineral (acc)         (19) mineral (acc) '
168      print*,'           (9) mineral (coa)         (20) mineral (coa) '
169      print*,'          (10) min.tra. (neu)        (21) min.tra. (neu)'
170      print*,'          (11) sulfate               (22) sulfate       '
171      print*,'    '
172      print*,'          (12) Aerosol Types         (23) Profil Typ'
173      write(*,111)
174  111 format(/,' please select number (0=END): ')
175
176      read(*,*) ii
177      if (ii.lt.0.or.ii.gt.23) then
178         print*,'       wrong value! Try Again!'
179         goto 1100
180      else if (ii.eq.0) then
181         ende=.true.
182         goto 1100
183      else
184      end if
185
186      parnam=parn(ii)
187
188      read (7,*) imonat
189      nx=0
190      ny=0
191      do iy=1,37
192         do ix=1,72
193            t4=0
194            rm4=0.
195            t5=0
196            rm5=0.
197            read (7,200,end=99) yi,xi,nl,np,at,rh,dn,nc,
198     *                          t1,rm1,t2,rm2,t3,rm3
199            if (nc.gt.3) read (7,202) t4,rm4
200            if (nc.gt.4) then
201               print*,' ACHTUNG: mehr als 4 Komponenten in 1. Schicht'
202               print*,' Breite:',yi,' L„nge:',xi
203               stop
204            end if
205            if (nl.eq.2) then
206               read (7,203) dn2,nc2,t5,rm5
207               if (nc2.gt.1) then
208                  print*,' ACHTUNG: ',nc2,' Komponenten in',
209     *                   ' zweiter Schicht bei ',yi,xi
210                  stop
211               end if
212            else if (nl.lt.1.or.nl.gt.2) then
213               print*,' ACHTUNG:',nl,' Schichten bei',yi,xi
214               stop
215            end if
216            if (ii.eq.1) then                 ! number density
217               w(ix,iy)=dn
218               if (nl.eq.2) then
219                  w2(ix,iy)=dn2
220               else
221                  w2(ix,iy)=0.
222               end if
223            else if (ii.eq.2) then            ! insoluble
224               if(t1.eq.1) then
225                  w(ix,iy)=rm1
226               else if(t2.eq.1) then
227                  w(ix,iy)=rm2
228               else if(t3.eq.1.and.nc.gt.2) then
229                  w(ix,iy)=rm3
230               else if(t4.eq.1.and.nc.gt.3) then
231                  w(ix,iy)=rm4
232               else
233                  w(ix,iy)=0.
234               end if
235               w(ix,iy)=w(ix,iy)*dn
236            else if (ii.eq.3) then            ! watersoluble
237               if(t1.eq.2) then
238                  w(ix,iy)=rm1
239               else if(t2.eq.2) then
240                  w(ix,iy)=rm2
241               else if(t3.eq.2.and.nc.gt.2) then
242                  w(ix,iy)=rm3
243               else if(t4.eq.2.and.nc.gt.3) then
244                  w(ix,iy)=rm4
245               else
246                  w(ix,iy)=0.
247               end if
248               w(ix,iy)=w(ix,iy)*dn
249            else if (ii.eq.4) then            ! soot
250               if(t1.eq.3) then
251                  w(ix,iy)=rm1
252               else if(t2.eq.3) then
253                  w(ix,iy)=rm2
254               else if(t3.eq.3.and.nc.gt.2) then
255                  w(ix,iy)=rm3
256               else if(t4.eq.3.and.nc.gt.3) then
257                  w(ix,iy)=rm4
258               else
259                  w(ix,iy)=0.
260               end if
261               w(ix,iy)=w(ix,iy)*dn
262            else if (ii.eq.5) then            ! sea salt (acc.)
263               if(t1.eq.4) then
264                  w(ix,iy)=rm1
265               else if(t2.eq.4) then
266                  w(ix,iy)=rm2
267               else if(t3.eq.4.and.nc.gt.2) then
268                  w(ix,iy)=rm3
269               else if(t4.eq.4.and.nc.gt.3) then
270                  w(ix,iy)=rm4
271               else
272                  w(ix,iy)=0.
273               end if
274               w(ix,iy)=w(ix,iy)*dn
275            else if (ii.eq.6) then            ! sea salt (coa.)
276               if(t1.eq.5) then
277                  w(ix,iy)=rm1
278               else if(t2.eq.5) then
279                  w(ix,iy)=rm2
280               else if(t3.eq.5.and.nc.gt.2) then
281                  w(ix,iy)=rm3
282               else if(t4.eq.5.and.nc.gt.3) then
283                  w(ix,iy)=rm4
284               else
285                  w(ix,iy)=0.
286               end if
287               w(ix,iy)=w(ix,iy)*dn
288            else if (ii.eq.7) then            ! mineral (nuc.)
289               if(t1.eq.6) then
290                  w(ix,iy)=rm1
291               else if(t2.eq.6) then
292                  w(ix,iy)=rm2
293               else if(t3.eq.6.and.nc.gt.2) then
294                  w(ix,iy)=rm3
295               else if(t4.eq.6.and.nc.gt.3) then
296                  w(ix,iy)=rm4
297               else
298                  w(ix,iy)=0.
299               end if
300               w(ix,iy)=w(ix,iy)*dn
301            else if (ii.eq.8) then            ! mineral (acc.)
302               if(t1.eq.7) then
303                  w(ix,iy)=rm1
304               else if(t2.eq.7) then
305                  w(ix,iy)=rm2
306               else if(t3.eq.7.and.nc.gt.2) then
307                  w(ix,iy)=rm3
308               else if(t4.eq.7.and.nc.gt.3) then
309                  w(ix,iy)=rm4
310               else
311                  w(ix,iy)=0.
312               end if
313               w(ix,iy)=w(ix,iy)*dn
314            else if (ii.eq.9) then            ! mineral (coa.)
315               if(t1.eq.8) then
316                  w(ix,iy)=rm1
317               else if(t2.eq.8) then
318                  w(ix,iy)=rm2
319               else if(t3.eq.8.and.nc.gt.2) then
320                  w(ix,iy)=rm3
321               else if(t4.eq.8.and.nc.gt.3) then
322                  w(ix,iy)=rm4
323               else
324                  w(ix,iy)=0.
325               end if
326               w(ix,iy)=w(ix,iy)*dn
327            else if (ii.eq.10) then            ! mineral (tra.)
328               if(t1.eq.9) then
329                  w(ix,iy)=rm1
330               else if(t2.eq.9) then
331                  w(ix,iy)=rm2
332               else if(t3.eq.9.and.nc.gt.2) then
333                  w(ix,iy)=rm3
334               else if(t4.eq.9.and.nc.gt.3) then
335                  w(ix,iy)=rm4
336               else if(t5.eq.9) then
337                  w(ix,iy)=rm5
338                  dn=dn2
339               else
340                  w(ix,iy)=0.
341               end if
342               w(ix,iy)=w(ix,iy)*dn
343            else if (ii.eq.11) then            ! sulfate
344               if(t1.eq.10) then
345                  w(ix,iy)=rm1
346               else if(t2.eq.10) then
347                  w(ix,iy)=rm2
348               else if(t3.eq.10.and.nc.gt.2) then
349                  w(ix,iy)=rm3
350               else if(t4.eq.10.and.nc.gt.3) then
351                  w(ix,iy)=rm4
352               else
353                  w(ix,iy)=0.
354               end if
355               w(ix,iy)=w(ix,iy)*dn
356            else if (ii.eq.12) then            ! Aerosol Types
357               if (at.eq.'CC'.or.at.eq.'RU') then
358                  w(ix,iy)=1.
359               else if (at.eq.'CA') then
360                  w(ix,iy)=2.
361               else if (at.eq.'MI') then
362                  w(ix,iy)=3.
363               else if (at.eq.'UR') then
364                  w(ix,iy)=4.
365               else if (at.eq.'MC'.and.nl.eq.1) then
366                  w(ix,iy)=5.
367               else if (at.eq.'MC'.and.nl.eq.2) then
368                  w(ix,iy)=6.
369               else if (at.eq.'MP'.and.nl.eq.2) then
370                  w(ix,iy)=7.
371               else if (at.eq.'MP'.and.nl.eq.1) then
372                  w(ix,iy)=8.
373               else if (at.eq.'NP') then
374                  w(ix,iy)=9.
375               else if (at.eq.'SP') then
376                  w(ix,iy)=10.
377               else
378                  print*,'Aerosoltyp ',at,' entdeckt'
379                  w(ix,iy)=11.
380               end if
381            else if (ii.eq.13) then            ! insoluble
382               if(t1.eq.1) then
383                  w(ix,iy)=rm1
384               else if(t2.eq.1) then
385                  w(ix,iy)=rm2
386               else if(t3.eq.1.and.nc.gt.2) then
387                  w(ix,iy)=rm3
388               else if(t4.eq.1.and.nc.gt.2) then
389                  w(ix,iy)=rm4
390               else
391                  w(ix,iy)=0.
392               end if
393               w(ix,iy)=w(ix,iy)*dn*vn(1)*dens(1)*10.**(-6)
394            else if (ii.eq.14) then            ! watersoluble
395               if(t1.eq.2) then
396                  w(ix,iy)=rm1
397               else if(t2.eq.2) then
398                  w(ix,iy)=rm2
399               else if(t3.eq.2.and.nc.gt.2) then
400                  w(ix,iy)=rm3
401               else if(t4.eq.2.and.nc.gt.2) then
402                  w(ix,iy)=rm4
403               else
404                  w(ix,iy)=0.
405               end if
406               w(ix,iy)=w(ix,iy)*dn*vn(2)*dens(2)*10.**(-6)
407            else if (ii.eq.15) then            ! soot
408               if(t1.eq.3) then
409                  w(ix,iy)=rm1
410               else if(t2.eq.3) then
411                  w(ix,iy)=rm2
412               else if(t3.eq.3.and.nc.gt.2) then
413                  w(ix,iy)=rm3
414               else if(t4.eq.3.and.nc.gt.2) then
415                  w(ix,iy)=rm4
416               else
417                  w(ix,iy)=0.
418               end if
419               w(ix,iy)=w(ix,iy)*dn*vn(3)*dens(3)*10.**(-6)
420            else if (ii.eq.16) then            ! sea salt (acc.)
421               if(t1.eq.4) then
422                  w(ix,iy)=rm1
423               else if(t2.eq.4) then
424                  w(ix,iy)=rm2
425               else if(t3.eq.4.and.nc.gt.2) then
426                  w(ix,iy)=rm3
427               else if(t4.eq.4.and.nc.gt.2) then
428                  w(ix,iy)=rm4
429               end if
430               w(ix,iy)=w(ix,iy)*dn*vn(4)*dens(4)*10.**(-6)
431            else if (ii.eq.17) then            ! sea salt (coa.)
432               if(t1.eq.5) then
433                  w(ix,iy)=rm1
434               else if(t2.eq.5) then
435                  w(ix,iy)=rm2
436               else if(t3.eq.5.and.nc.gt.2) then
437                  w(ix,iy)=rm3
438               else if(t4.eq.5.and.nc.gt.2) then
439                  w(ix,iy)=rm4
440               else
441                  w(ix,iy)=0.
442               end if
443               w(ix,iy)=w(ix,iy)*dn*vn(5)*dens(5)*10.**(-6)
444            else if (ii.eq.18) then            ! mineral (nuc.)
445               if(t1.eq.6) then
446                  w(ix,iy)=rm1
447               else if(t2.eq.6) then
448                  w(ix,iy)=rm2
449               else if(t3.eq.6.and.nc.gt.2) then
450                  w(ix,iy)=rm3
451               else if(t4.eq.6.and.nc.gt.2) then
452                  w(ix,iy)=rm4
453               else
454                  w(ix,iy)=0.
455               end if
456               w(ix,iy)=w(ix,iy)*dn*vn(6)*dens(6)*10.**(-6)
457            else if (ii.eq.19) then            ! mineral (acc.)
458               if(t1.eq.7) then
459                  w(ix,iy)=rm1
460               else if(t2.eq.7) then
461                  w(ix,iy)=rm2
462               else if(t3.eq.7.and.nc.gt.2) then
463                  w(ix,iy)=rm3
464               else if(t4.eq.7.and.nc.gt.2) then
465                  w(ix,iy)=rm4
466               else
467                  w(ix,iy)=0.
468               end if
469               w(ix,iy)=w(ix,iy)*dn*vn(7)*dens(7)*10.**(-6)
470            else if (ii.eq.20) then            ! mineral (coa.)
471               if(t1.eq.8) then
472                  w(ix,iy)=rm1
473               else if(t2.eq.8) then
474                  w(ix,iy)=rm2
475               else if(t3.eq.8.and.nc.gt.2) then
476                  w(ix,iy)=rm3
477               else if(t4.eq.8.and.nc.gt.2) then
478                  w(ix,iy)=rm4
479               else
480                  w(ix,iy)=0.
481               end if
482               w(ix,iy)=w(ix,iy)*dn*vn(8)*dens(8)*10.**(-6)
483            else if (ii.eq.21) then            ! mineral (tra.)
484               if(t1.eq.9) then
485                  w(ix,iy)=rm1
486               else if(t2.eq.9) then
487                  w(ix,iy)=rm2
488               else if(t3.eq.9.and.nc.gt.2) then
489                  w(ix,iy)=rm3
490               else if(t4.eq.9.and.nc.gt.2) then
491                  w(ix,iy)=rm4
492               else if(t5.eq.9) then
493                  w(ix,iy)=rm5
494                  dn=dn2
495               else
496                  w(ix,iy)=0.
497               end if
498               w(ix,iy)=w(ix,iy)*dn*vn(9)*dens(9)*10.**(-6)
499            else if (ii.eq.22) then            ! sulfate
500               if(t1.eq.10) then
501                  w(ix,iy)=rm1
502               else if(t2.eq.10) then
503                  w(ix,iy)=rm2
504               else if(t3.eq.10.and.nc.gt.2) then
505                  w(ix,iy)=rm3
506               else if(t4.eq.10.and.nc.gt.2) then
507                  w(ix,iy)=rm4
508               else
509                  w(ix,iy)=0.
510               end if
511               w(ix,iy)=w(ix,iy)*dn*vn(10)*dens(10)*10.**(-6)
512            else if (ii.eq.23) then            ! Profil Typ
513               w(ix,iy)=np
514            end if
515
516            if (iy.eq.1) then
517               nx=nx+1
518               x(ix)=xi
519            end if
520            if (ix.eq.1) then
521               ny=ny+1
522               y(iy)=yi
523            end if
524         end do
525      end do
526   99 continue
527
528ccccc -----------------------------------------------------------------c
529c     output                                                           c
530ccccc -----------------------------------------------------------------c
531
532      file=parout(ii)//'.out'
533      open(9,file=file)
534     
535      write (9,160) season,parnam,parout(ii)
536
537      do iy=1,37
538         do ix=1,72
539            write (9,150) y(iy),x(ix),w(ix,iy)
540         end do
541      end do
542 
543      close (9)
544
545      rewind (7)
546
547      end do      ! Ende der Endlosschleife vom Anfang
548
549      close (7)
550
551ccccc -----------------------------------------------------------------c
552c     Formate                                                          c
553ccccc -----------------------------------------------------------------c
554
555  100 format (8e10.3)
556  101 format (i2)
557  102 format ('1')
558  103 format (' '/' '/' '/' ')
559  104 format (' '/' ')
560 
561  150 format (1x,2i5,1pe10.3)
562  160 format ('# Global Aerosol Data Set, Version 2.2a',/
563     *        '#',/
564     *        '# ',a7,/
565     *        '# ',/
566     *        '# ',a31,/,
567     *        '# ',/
568     *        '#  LAT  LON ',a5)
569     
570  200 format (2i4,2i3,1x,a2,4x,i2,e10.3,3x,i1,3(i3,e10.3))
571  202 format (37x,3(i3,e10.3))
572  203 format (23x,e10.3,3x,i1,1(i3,e10.3))
573  222 format (i5,1x,13e9.3,i5)
574  223 format (6x,13(3x,i3,3x))
575  224 format (6x,37(i3))
576  225 format (i5,1x,13(3x,a3,3x),i5)
577  227 format (i5,1x,11(3x,a3,3x),i5)
578  228 format (i5,1x,11e9.3,i5)
579  229 format (6x,13e9.3)
580  249 format ('1',' TAPE201.RDT am  ',i2,'.',i2,'.',i4,2x,
581     *        i2,':',i2,':',i2)
582  250 format ('1',' TAPE207.RDT am  ',i2,'.',i2,'.',i4,2x,
583     *        i2,':',i2,':',i2)
584  239 format (' TAPE201.RDT am  ',i2,'.',i2,'.',i4,2x,
585     *        i2,':',i2,':',i2)
586  240 format (' TAPE207.RDT am  ',i2,'.',i2,'.',i4,2x,
587     *        i2,':',i2,':',i2)
588  251 format (1x,a31)
589  300 format (a1)
590  333 format (i2,'.',i2,'.',i4,2x,i2,':',i2,':',i2)
591  350 format (60x,a10,2x,a8)
592  360 format (a10,2x,a8)
593  370 format ('Minimum: ',1pe10.3,'   Maximum: ',e10.3)
594  400 format (1x,10i2)
595  500 format (1pe8.2)
596  600 format (9x,a7)
597  650 format (9x,a50)
598  660 format (25x,i4,22x,i4,17x,i4)
599  700 format (71x,i1)
600  800 format (a7)
601  900 format (13x,f6.3,22x,a2)
602  950 format (20x,'  relative humidity: ',a2,' %')
603  960 format (a31)
604  980 format ('GADS')
605
606      stop
607      end
608      subroutine opt
609
610ccccc -----------------------------------------------------------------c
611c      Berechnung der optischen Aerosoldaten aus den mikrophysika-     c
612c      lischen Rohdaten.                                               c
613c                                                                      c
614c      ATLOPT ist eine modifizierte Version des Programms OPAC und     c
615c      steuert den Aufruf der Unterprogramme:                          c
616c                                                                      c
617c      - HEAD4                                                         c
618c      - LOCATE                                                        c
619c      - D4RAW                                                         c
620c      - OPTCOM                                                        c
621c      - OPTPAR                                                        c
622c      - OUT4                                                          c
623c                                                                      c
624c      Diese Unterprogramme sind an ATLOPT angeh„ngt.                  c
625c                                                                      c
626c      ab 14.12.92 anderes Format der neuen Mie-Rechnungen eingebaut   c
627c      ab 04.11.93 neue Parameter scattering, absorption, omega ratio  c
628c      ab 13.05.94 OPTCOM: Russ quillt nicht mehr                      c
629c      ab 27.07.94 opt. Dicke fuer alle Wellenlaengen berechnet        c
630c                                                                      c
631c      18.11.97 GADS 2.1                                               c
632c                                                                      c
633c      18.11.97                                                M. Hess c       
634ccccc -----------------------------------------------------------------c
635
636      integer   prnr,acnr,njc,rht
637      real      n,numden
638
639      character*1 ws,dum
640      character*2 chum
641      character*3  atn,pat
642      character*8  opanam,optnam
643      character*4  comnam
644      character*7  cseas
645      character*11 tseas
646      character*20 catyp
647      character*30 typnam
648      character*50 area
649     
650      common /prog/   nprog
651      common /profi/  nil(10),hfta(10),hstra(10),
652     *                h0(2,10),h1(2,10),hm(2,10)
653      common /buffer/ ibuf,kbuf(20),extbuf(20),scabuf(20),absbuf(20),
654     *                sisbuf(20),asybuf(20),bacbuf(20),phabuf(112,20),
655     *                brebuf(20),bimbuf(20),mbuf
656      common /mipoi/  latx,lonx,nl,prnr,rht(2),n(2),
657     *                njc(2),acnr(5,2),acmr(5,2),nh(2),atn(2),pat(2)
658      common /oppoi/  ext(2,5),sca(2,5),abs(2,5),sis(2,5),asy(2,5),
659     *                bac(2,5),pha(2,5,112),bre(2,5),bim(2,5)
660
661      common /atyp/   natyp,mcomp,ncomp(10),numden,
662     *                catyp,typnam,comnam(20)
663      common /layer/  nlay,mlay,nltyp(10),parlay(10,2),boundl(10),
664     *                boundu(10)
665      common /param/  nptyp,mpar,par(5)
666      common /season/ nseas,cseas,tseas
667      common /opar/   mopar,jnopar(13),nop,opanam(13),optnam(13)
668      common /wavel/  mlamb,alamb(61),niw
669      common /hum/    khum(8),ahum(8),nih,nhum(8),mhum,chum(8)
670      common /geog/   lata,late,lati,lona,lone,loni,na,area
671      common /norm/   norm,mixnor
672       
673      data alamb /0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,
674     *            0.9,1.0,1.25,1.5,1.75,2.0,2.5,3.0,3.2,3.39,3.5,3.75,
675     *            4.0,4.5,5.0,5.5,6.0,6.2,6.5,7.2,7.9,8.2,8.5,8.7,9.0,
676     *            9.2,9.5,9.8,10.0,10.6,11.0,11.5,12.5,13.0,14.0,14.8,
677     *            15.0,16.4,17.2,18.0,18.5,20.0,21.3,22.5,25.0,27.9,30.,
678     *            35.0,40.0/,mlamb/61/
679     
680      data ahum /0.,50.,70.,80.,90.,95.,98.,99./
681      data nhum /0,50,70,80,90,95,98,99/,mhum/8/
682      data chum /'00','50','70','80','90','95','98','99'/
683
684      data comnam /'inso','waso','soot','ssam','sscm','minm','miam',
685     *             'micm','mitr','suso','stco','stma','cucc','cucp',
686     *             'cuma','fog-','cir1','cir2','cir3','    '/
687
688      data optnam /'ext.coef','sca.coef','abs.coef','sisc.alb',
689     *             'asym.par','op.depth',
690     *             '        ','turb.fac','li.ratio','pha.func',
691     *             'ext.rat ','abs.rat ',
692     *             '        '/
693     
694      data jnopar/1,1,1,1,1,1,0,0,1,0,0,0,0/,nop/7/
695
696CCCCC -----------------------------------------------------------------C
697c     some definitions for this version                                c
698CCCCC -----------------------------------------------------------------C
699       
700      niw=1
701      njh=1
702      nih=1
703      lata=90
704      late=-90
705      lati=5
706      lona=-180
707      lone=175
708      loni=5
709      nlmal=1
710     
711      norm=1
712      nprog=4
713
714      ntape=22
715     
716      ip=0
717      do i=1,13
718         if (jnopar(i).eq.1) then
719            ip=ip+1
720            opanam(ip)=optnam(i)
721         end if
722      end do
723
724ccccc -----------------------------------------------------------------c
725c     Abfrage, was geplottet werden soll                               c
726ccccc -----------------------------------------------------------------c
727
728 1001 print*,'  '
729      write(*,154)
730  154 format(' (w)inter or (s)ummer? ')
731      read (*,'(a)') ws
732      if (ws.eq.'w') then
733         open(ntape,file='../glodat/winter.dat')
734         read (ntape,'(a1)') dum
735         cseas='winter '
736      else if (ws.eq.'s') then
737         open(ntape,file='../glodat/summer.dat')
738         read (ntape,'(a1)') dum
739         cseas='summer '
740      else
741         print*,' wrong input! try again!'
742         goto 1001
743      end if
744
745ccccc ----------------------------------------------------------------c
746c     Input: wavelength                                               c
747ccccc ----------------------------------------------------------------c
748
749      print*,' please select wavelength: '
750      print*,' '
751     
752      nwel=mlamb
753      do 11 iwel=1,22
754         if (nwel.ge.(iwel+44)) then
755            write(*,114) iwel,alamb(iwel),(iwel+22),alamb(iwel+22),
756     *                   (iwel+44),alamb(iwel+44)
757         else if (nwel.ge.(iwel+22)) then
758            write(*,113) iwel,alamb(iwel),(iwel+22),alamb(iwel+22)
759         else
760            write(*,111) iwel,alamb(iwel)
761         end if
762   11 continue
763
764  111 format(5x,'(',i2,')',3x,f5.2,1x,'um')
765  113 format(5x,'(',i2,')',3x,f5.2,1x,'um',5x,'(',i2,')',
766     *       3x,f5.2,1x,'um')
767  114 format(5x,'(',i2,')',3x,f5.2,1x,'um',5x,'(',i2,')',
768     *       3x,f5.2,1x,'um',
769     *       5x,'(',i2,')',3x,f5.2,1x,'um')
770
771  909 write (*,*) '?'
772      read (*,*) iwel
773      if (iwel.lt.1.or.iwel.gt.nwel) then
774         print*,' wrong number! try again! '
775         goto 909
776      end if
777
778ccccc ----------------------------------------------------------------c
779c     Input: humidity                                                 c
780ccccc ----------------------------------------------------------------c
781
782      print*,' please select rel. humidity: '
783      print*,' '
784      do ihum=1,mhum
785         write(*,121) ihum,nhum(ihum)
786  121    format (5x,'(',i2,')',3x,i2,' %')       
787      end do
788     
789  908 write (*,*) '?'
790      read (*,*) ihum
791      if (ihum.lt.1.or.ihum.gt.mhum) then
792         print*,' wrong number! try again! '
793         goto 908
794      end if
795
796CCCCC -----------------------------------------------------------------C
797C     EINLESEN DER HOEHEN-PROFILE vom File TAPE9                       c
798C                                                                      C
799C     HM    : EFFEKTIVE SCHICHTDICKE (HOMOGENE VERTEILUNG)             C
800C     HFTA  : SCHICHTDICKE DES FREIEN TROP. AEROSOLS IN KM             C
801C     HSTRA : SCHICHTDICKE DES STRATOSPH. AEROSOLS IN KM               C
802C     NIL   : ANZAHL DER SCHICHTEN                                     C
803CCCCC -----------------------------------------------------------------C
804
805c      print*,' Anfang prof'
806       call prof
807c      print*,' Ende prof'
808
809ccccc -----------------------------------------------------------------c
810c      Beschriftung des Output-Files                                   c
811ccccc -----------------------------------------------------------------c
812
813c      print*,' Anfang head4'
814       call head4 (iwel,ihum)
815c      print*,' Ende head4'
816
817ccccc -----------------------------------------------------------------c
818c      Schleife ber alle verlangten Wellenl„ngen und Feuchteklassen   c
819ccccc -----------------------------------------------------------------c
820
821       do il=1,niw
822
823          do ih=1,njh
824
825ccccc -----------------------------------------------------------------c
826c      Schleife ber alle verlangten geographischen Koordinaten        c
827ccccc -----------------------------------------------------------------c
828
829             do ilat=lata,late,-lati
830             
831                do ilmal=1,nlmal
832                   do ilon=lona,lone,loni
833                   
834ccccc -----------------------------------------------------------------c
835c      Einlesen der Rohdaten von den Files TAPE201, TAPE207:           c
836c     ------------------------------------------------------           c
837C     LAT       : LATITUDE                                             C
838C     LON       : LONGITUDE                                            C
839C     NL        : NUMBER OF AEROSOL LAYERS                             C
840C                 (=2 FOR MARITIME-MINERAL,=1 FOR MARI.)               C
841C     PRNR      : PROFIL NUMBER                                        C
842C     NT        : NUMBER OF AEROSOL TYPE                               C
843C     PAT       : AEROSOL PROFIL TYPE                                  C
844C     NH        : NUMBER OF REL. HUMIDITY CLASS                        C
845C     N         : TOTAL NUMBER CONCENTRATION                           C
846C     NJC       : NUMBER OF AEROSOL COMPONENT                          C
847C     ACNR      : AEROSOL COMPONENT                                    C
848C     ACMR      : MIXING RATIO                                         C
849C                 (PARTIAL NUMBER CONCENTRATION/TOTAL NUMBER CONC.)    C
850ccccc -----------------------------------------------------------------c
851
852c       print*,' Anfang d4raw'
853                      call d4raw (ilat,ilon,ntape)
854c       print*,' Ende d4raw'
855
856ccccc -----------------------------------------------------------------c
857c      Einlesen der optischen Rohdaten von den Files winter.dat and    c
858c      summer.dat                                                      c
859ccccc -----------------------------------------------------------------c
860
861c       print*,' Anfang optcom'
862                      call optcom (iwel,ihum)
863c       print*,' ende optcom'
864
865ccccc -----------------------------------------------------------------c
866c      Berechnung der optischen Parameter am aktuellen Gitterpunkt     c
867ccccc -----------------------------------------------------------------c
868
869c       print*,' Anfang optpar',ilat,ilon
870                      call optpar(iwel,ihum)
871c       print*,' Ende optpar',ilat,ilon
872
873                   end do
874                end do
875             end do
876          end do
877       end do
878       
879       close (ntape)
880       close (10)
881
882       return
883       end
884
885CCCCC *****************************************************************C
886      SUBROUTINE PROF
887C     *****************************************************************C
888C                                                                      C
889C     -----------------------------------------------------------------C
890C     EINLESEN DER HOEHEN-PROFILE vom File profiles.dat und der        C
891C     Extinktionskoeffizienten der oberen Atmosph„re von extcof.dat    C
892CCCCC -----------------------------------------------------------------C
893
894      common /wavel/  mlamb,alamb(61),niw
895      common /profi/ nil(10),hfta(10),hstra(10),
896     *                h0(2,10),h1(2,10),hm(2,10)
897      COMMON /FTASTR/ EXTFTA(61),EXTSTR(61)
898
899CCCCC -----------------------------------------------------------------C
900C     ES GIBT 7 PROFIL-TYPEN. Folgende Daten werden eingelesen:        C
901c                                                                      c
902c           iip: Nummer des Profiltyps                                 c
903c       nil(ip): Zahl der Schichten fuer Typ ip (wie in tape201 usw.)  c
904c      hfta(ip): Hoehe der Schicht fuer das freie troposph. Aerosol    c
905c     hstra(ip): Hoehe der Schicht des stratosphaerischen Aerosols     c
906c     h0(il,ip): Untergrenze der Schicht il                            c
907c     h1(il,ip): Obergrenze der Schicht                                c
908c     hm(il,ip): effektive Dicke der Schicht il fuer den Typ ip        c
909CCCCC -----------------------------------------------------------------C
910
911      open (8,file='profiles.dat')
912
913      nprof=7
914      DO IP=1,nprof
915         READ(8,8010) IIP,NIL(IP),HFTA(IP),HSTRA(IP)
916         READ(8,8020) (H0(IL,IP),H1(IL,IP), HM(IL,IP),IL=1,NIL(IP) )
917 8010    FORMAT(I3,I3,2F8.2)
918 8020    FORMAT(2F5.1,F10.3)
919      end do
920
921      close (8)
922
923CCCCC -----------------------------------------------------------------C
924C     Einlesen der Extinktionskoeffizienten fr die obere Atmosph„re:  C
925C                                                                      C
926C     EXTINCTION COEFFICIENT -  FREE TROPOSPHERIC AEROSOL  +           C
927C     EXTINCTION COEFFICIENT -    STRATOSPHERIC   AEROSOL              C
928C                                                                      C
929C     UEBERSPRINGEN DER ERSTEN BEIDEN ZEILEN von TAPE9                 C
930CCCCC -----------------------------------------------------------------C
931
932      open (9,file='extback.dat')
933      IL=1
934      READ(9,'(/)')
935      DO IWL=1,mlamb
936         READ(9,*) WAVE,EXTFT,EXTST
937c        do ila=1,niw
938c           IF (WAVE.EQ.alamb(ila)) THEN
939               EXTFTA(IWL)=EXTFT
940               EXTSTR(IWL)=EXTST
941c              IL=IL+1
942c           END IF
943c        end do
944      end do
945
946      close (9)
947
948      RETURN
949      END
950
951ccccc *****************************************************************c
952      subroutine head4 (il,ih)
953c     *****************************************************************c
954c                                                                      c
955c     -----------------------------------------------------------------c
956c     Beschriftung des Output-Files TAPE10                                 c
957ccccc -----------------------------------------------------------------c
958
959      real mixrat,numden
960     
961      character*2 chum
962      character*4 comnam
963      character opanam*8,cseas*7,tseas*11,optnam*8,opnam(10)*8
964      character catyp*20,area*50,typnam*30
965
966      common /season/ nseas,cseas,tseas
967      common /geog/   lata,late,lati,lona,lone,loni,na,area
968      common /hum/    khum(8),ahum(8),nih,nhum(8),mhum,chum(8)
969      common /norm/   norm,mixnor
970      common /numdis/ sigma(10),rmin(10,8),rmax(10,8),rmod(10,8),
971     *                mixrat(10),dens(10,8)
972      common /atyp/   natyp,mcomp,ncomp(10),numden,
973     *                catyp,typnam,comnam(20)
974      common /opar/   mopar,jnopar(13),nop,opanam(13),optnam(13)
975      common /wavel/  mlamb,alamb(61),niw
976      common /angle/  jnangle(112),angle(112),nia
977
978CCCCC -----------------------------------------------------------------C
979C     Kopf des output-files                                            C
980CCCCC -----------------------------------------------------------------C
981
982c      if (ih.eq.1.and.il.eq.1) then
983         open (10,file='aererg')
984
985         write(10,100) cseas
986  100      format('# Global Aerosol Data Set, Version 2.2a'/,
987     *          '#'/
988     *          '# ',a7,/
989     *          '#')
990  107    format('====================================================',
991     *          '============================')
992c      end if
993
994CCCCC -----------------------------------------------------------------C
995C     Beschriftung fr verschiedene Wellenlaengen und rel. Feuchten    c               c
996CCCCC -----------------------------------------------------------------C
997
998      if (nih.ne.0) then
999         WRITE(10,4000) alamb(il),ahum(ih)
1000 4000    FORMAT('#'/
1001     *          '# wavelength: ',f6.3,3x,'relative humidity: ',F3.0,'%'/
1002     *          '#')
1003      else
1004         WRITE(10,4004) alamb(il)
1005 4004    FORMAT(/' wavelength: ',f6.3,3x,'relative humidity: -- %')
1006      end if
1007
1008      if (jnopar(10).eq.1) then
1009         kop=nop-1
1010      else
1011         kop=nop
1012      end if
1013
1014      do iop=1,kop
1015         opnam(iop)=opanam(iop)
1016      end do
1017
1018      if (kop.le.10) then
1019         write(10,4001) (opnam(in),in=1,kop)
1020 4001    format('# LAT  LON NL ',10(1x,a8,1x))
1021      else
1022         write(10,4001) (opnam(in),in=1,10)
1023         write(10,4002) (opnam(in),in=11,kop)
1024 4002    format('                               ',5(1x,a8,1x))
1025      end if
1026     
1027      write(10,4003)
1028 4003 format('#',13x,'  [1/km]  ','  [1/km]  ','  [1/km]  ',
1029     *       30x,'   [sr]')
1030
1031      return
1032      end
1033
1034CCCCC *****************************************************************C
1035      subroutine d4raw (lat,lon,ntape)
1036C     *****************************************************************C
1037C                                                                      C
1038C     -----------------------------------------------------------------C
1039C     EINLESEN DER DATEN VON DEN ROHDATEN-FILES TAPE201-TAPE212        C
1040CCCCC -----------------------------------------------------------------C
1041
1042      IMPLICIT CHARACTER*3 (Z)
1043
1044      integer prnr,acnr,njc,rht
1045      real n
1046      character*3 atn,pat
1047
1048      common /prog/  nprog
1049      common /mipoi/ latx,lonx,nl,prnr,rht(2),n(2),
1050     *               njc(2),acnr(5,2),acmr(5,2),nh(2),atn(2),pat(2)
1051      common /test/  itest(2)
1052      common /dat/   zat(11),zrh(8)
1053
1054  111 READ(NTAPE,1020,end=999) LATX,LONX,NL,PRNR,
1055     +   ATN(1),PAT(1),RHT(1),N(1),NJC(1),
1056     +      ( ACNR(JC,1),ACMR(JC,1),JC=1,3 )
1057 1020 FORMAT(2I4,2I3,(2A3,I3,E10.3,I4,3(I3,E10.4)))
1058
1059c      write(*,1020) LATX,LONX,NL,PRNR,
1060c     +   ATN(1),PAT(1),RHT(1),N(1),NJC(1),
1061c     +      ( ACNR(JC,1),ACMR(JC,1),JC=1,3 )
1062
1063      IF(LATX.NE.LAT .OR.  LONX.NE.LON) THEN
1064         print*,' Achtung, falsche Koordinaten: ',latx,lonx
1065         GOTO 111
1066      END IF
1067
1068      IF (NJC(1).GT.3) THEN
1069          READ(NTAPE,1025,end=999) ( ACNR(JC,1),ACMR(JC,1),JC=4,NJC(1))
1070 1025     FORMAT(37X,3(I3,E10.4))
1071c          write(*,1025) ( ACNR(JC,1),ACMR(JC,1),JC=4,NJC(1))
1072      END IF
1073
1074      IF(NL.NE.1) THEN
1075         DO 10 L=2,NL
1076         READ(NTAPE,1021,end=999)  ATN(L),PAT(L),RHT(L),N(L),NJC(L),
1077     +                   ( ACNR(JC,L),ACMR(JC,L),JC=1,3 )
1078 1021    FORMAT(14X,2A3,I3,E10.3,I4,5(I3,E10.4))
1079
1080         IF (NJC(L).GT.3) THEN
1081           READ(NTAPE,1025,end=999) ( ACNR(JC,L),ACMR(JC,L),JC=4,NJC(L))
1082         END IF
1083
1084   10    CONTINUE
1085      END IF
1086
1087      do l=1,nl
1088         sum=0.
1089         do ic=1,njc(l)
1090            sum=sum+acmr(ic,l)
1091         end do
1092         if (abs(sum-1.).ge.0.01) then
1093            print*,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
1094            print*,'***!!!    sum of mixing ratios is not 1.     !!!***'
1095            print*,'***!!! please have a look at errorfile *.err !!!***'
1096            print*,'***************************************************'
1097            write (2,1001) latx,lonx,sum
1098         end if
1099      end do
11001001  format (2i4,3x,1pe10.3)
1101
1102CCCCC -----------------------------------------------------------------C
1103C     BESTIMMUNG DER NUMMER DES AEROSOLTYPS UND DER FEUCHTEKLASSE      C
1104C     !!! AEROSOLTYP IST SCHICHTABHAENGIG (NOCH NICHT BERUECKSICHTIGT) C
1105CCCCC -----------------------------------------------------------------C
1106
1107c     DO 50 IT=1,11
1108c     IF(ATN(1).EQ.ZAT(IT)) THEN
1109c        NT=IT
1110c     END IF
1111c  50 CONTINUE
1112
1113      DO 60 IL=1,NL
1114      IF(RHT(IL).LE.30) THEN
1115         NH(IL)=1
1116      ELSE IF(RHT(IL).GT.30.AND.RHT(IL).LE.65) THEN
1117         NH(IL)=2
1118      ELSE IF(RHT(IL).GT.65.AND.RHT(IL).LE.75) THEN
1119         NH(IL)=3
1120      ELSE IF(RHT(IL).GT.75.AND.RHT(IL).LE.85) THEN
1121         NH(IL)=4
1122      ELSE IF(RHT(IL).GT.85.AND.RHT(IL).LE.92) THEN
1123         NH(IL)=5
1124      ELSE IF(RHT(IL).GT.92.AND.RHT(IL).LE.97) THEN
1125         NH(IL)=6
1126      ELSE IF(RHT(IL).EQ.98) THEN
1127         NH(IL)=7
1128      ELSE IF(RHT(IL).EQ.99) THEN
1129         NH(IL)=8
1130      END IF
1131   60 CONTINUE
1132
1133  999 RETURN
1134      END
1135
1136CCCCC *****************************************************************C
1137       subroutine optcom (ilamb,ihum)
1138C     *****************************************************************C
1139C                                                                      C
1140C     -----------------------------------------------------------------C
1141C      Einlesen der optischen Rohdaten in einen Puffer fr             C
1142C      20 Komponenten                                                  C
1143c                                                                      c
1144c     Bei den urspruenglichen Mie-Rechnungen sind alle Koeffizienten   c
1145c     und die Phasenfunktion in [1/cm] angegeben. Die neuen Rechnungen c
1146c     geben die Ergebnisse in [1/m]. Daher muessen die neuen           c
1147c     Ergebnisse in der ersten Zeile durch den Zusatz 'neu' gekenn-    c
1148c     zeichnet werden: z.B. TAPE741, neu                               c
1149c                                                                      c
1150c     13.05.94 Quellung von Russ ausgeschlossen.                       c
1151c     17.11.97 new file format in ../optdat/                           c 
1152c                                                                      c
1153c     Stand: 17.11.97                                         M. Hess  c
1154CCCCC -----------------------------------------------------------------C
1155
1156      integer prnr,acnr,njc,rht
1157      real n,numden
1158     
1159      character*1  dum
1160      character*2  chum
1161      character*3  atn,pat
1162      character*4  comnam
1163      character*8  opanam,optnam
1164      character*10 dum2
1165      character*18 tap
1166      character*20 catyp
1167      character*30 typnam
1168     
1169      logical exists,ende
1170
1171      common /opar/   mopar,jnopar(13),nop,opanam(13),optnam(13)
1172      common /wavel/  mlamb,alamb(61),niw
1173      common /hum/    khum(8),ahum(8),nih,nhum(8),mhum,chum(8)
1174      common /atyp/   natyp,mcomp,ncomp(10),numden,
1175     *                catyp,typnam,comnam(20)
1176      common /mipoi/  latx,lonx,nl,prnr,rht(2),n(2),
1177     *                njc(2),acnr(5,2),acmr(5,2),nh(2),atn(2),pat(2)
1178      common /oppoi/  ext(2,5),sca(2,5),abs(2,5),sis(2,5),asy(2,5),
1179     *              bac(2,5),pha(2,5,112),bre(2,5),bim(2,5)
1180      common /buffer/ ibuf,kbuf(20),extbuf(20),scabuf(20),absbuf(20),
1181     *                sisbuf(20),asybuf(20),bacbuf(20),phabuf(112,20),
1182     *            brebuf(20),bimbuf(20),mbuf
1183     
1184ccccc -----------------------------------------------------------------c
1185c      Schleife ber alle am Gitterpunkt vorkommenden Komponenten      c
1186ccccc -----------------------------------------------------------------c
1187
1188       do il=1,nl
1189       
1190          if (nih.eq.0) then
1191             do ihu=1,mhum
1192                if (nh(il).eq.ihu) then
1193                   khum(ihum)=ihu
1194                end if
1195             end do
1196          end if
1197
1198          do ic=1,njc(il)
1199         
1200c         print*,'Anfang Komponenten schleife: ',ic,njc(il)
1201         
1202             jc=acnr(ic,il)
1203             
1204ccccc -----------------------------------------------------------------c
1205c     Ausschluá der Quellung bei insoluble, Russ und den               c
1206c     mineralischen Komponenten und bei den Wolken                     c                          c
1207ccccc -----------------------------------------------------------------c
1208
1209            if ( jc.eq.1.or.jc.eq.3.or.(jc.ge.6.and.jc.le.9).or.
1210     *           jc.gt.10 ) then
1211               iht=1
1212               nta=700+(jc*10)+1
1213            else
1214c               iht=khum(ihum)
1215                iht=ihum
1216                nta=700+(jc*10)+iht
1217            end if
1218
1219ccccc -----------------------------------------------------------------c
1220c     Bestimmung des Filenamens der gesuchten Komponente aus           c
1221c     Komponentennummer und Feuchteklasse                              c
1222ccccc -----------------------------------------------------------------c
1223
1224            tap(1:10)='../optdat/'
1225            tap(11:14)=comnam(jc)
1226            tap(15:16)=chum(iht)
1227
1228            ntap=20
1229
1230ccccc -----------------------------------------------------------------c
1231c     Bestimmung der Kennnummer fr den Puffer ber die Wellenl„nge    c
1232c                                                                      c
1233c         nbuf: Kennummer der aktuellen Komponente fuer den Puffer     c
1234c     kbuf(20): Kennummern der gespeicherten Komponenten               c
1235c         mbuf: Position der aktuellen Komponente im Puffer            c
1236c         ibuf: Position bis zu der der Puffer belegt ist              c
1237ccccc -----------------------------------------------------------------c
1238
1239             nbuf=ilamb*1000+nta
1240c            print*,'nbuf= ',nbuf
1241
1242ccccc -----------------------------------------------------------------c
1243c      šberprfung des Puffers auf šbereinstimmung mit nbuf            c
1244ccccc -----------------------------------------------------------------c
1245
1246             exists=.false.
1247             do ib=1,20
1248                if (nbuf.eq.kbuf(ib)) then
1249                   exists=.true.
1250                   mbuf=ib
1251                   goto 10
1252                end if
1253             end do
1254   10      continue
1255c           print*,'mbuf= ',mbuf
1256
1257ccccc -----------------------------------------------------------------c
1258c      Einlesen der Komponentendaten, falls sie nicht im Puffer        c
1259c      stehen, sonst šbernahme aus dem Puffer                          c
1260ccccc -----------------------------------------------------------------c
1261
1262c       print*,exists
1263
1264         if (exists) then
1265            ext(il,ic)=extbuf(mbuf)
1266            sca(il,ic)=scabuf(mbuf)
1267            abs(il,ic)=absbuf(mbuf)
1268            sis(il,ic)=sisbuf(mbuf)
1269            asy(il,ic)=asybuf(mbuf)
1270            bac(il,ic)=bacbuf(mbuf)
1271            bre(il,ic)=brebuf(mbuf)
1272            bim(il,ic)=bimbuf(mbuf)
1273         else
1274            if (ibuf.lt.20) then
1275               ibuf=ibuf+1
1276            else
1277               print*,' ibuf= ',ibuf
1278               ibuf=1
1279            end if
1280c            print*,ibuf
1281           
1282            kbuf(ibuf)=nbuf
1283            open (ntap,file=tap,iostat=ios)
1284c            print*,'opened file ',tap,iostat
1285           
1286            if (ios.ne.0) then
1287               print*,' error while opening file ',tap
1288               print*,'  latitude: ',latx
1289               print*,' longitude: ',lonx
1290               stop
1291            end if
1292           
1293           
1294c ALTER INPUT           
1295c            read (ntap,200) dum
1296c            read (ntap,'(22(/))')
1297c            rlamb=0.
1298c            do while (rlamb.ne.alamb(ilamb))
1299c            do ila=1,ilamb
1300c               read (ntap,500) rlamb,extco,sisca,asymf,exn,refr,refi
1301c  500          format(2x,8e10.3)
1302c            end do
1303c            do ila=ilamb+1,mlamb
1304c               read (ntap,500) rl
1305c               print*,rl
1306c            end do
1307c                       
1308c            read (ntap,'(4(/))')
1309c           
1310c            ntheta=96
1311c            do it=1,ntheta
1312c           
1313c               read (ntap,510,end=511)
1314c     *               thet,(pha(il,ic,it),ila=1,ilamb)
1315c  510          format(1x,70e10.3)
1316c 
1317c               print*,it,thet,pha(il,ic,it)
1318c 
1319c            end do
1320c  511       continue
1321c 
1322c   
1323c  ENDE ALTER INPUT
1324 
1325         do iline=1,100
1326            read (ntap,220) dum2
1327            if (dum2.eq.'# optical ') then
1328               goto 2002
1329            end if
1330         end do
1331 2002    continue
1332         do iline=1,5
1333            read (ntap,200) dum
1334         end do                 
1335               
1336         do ila=1,mlamb
1337            read (ntap,500) rlamb,extco,scaco,absco,sisca,asymf,
1338     *                      exn,refr,refi
1339  500       format(2x,7e10.3,2e11.3)
1340 
1341            if (rlamb.eq.alamb(ilamb)) then
1342               ext(il,ic)=extco
1343               sca(il,ic)=scaco
1344               abs(il,ic)=absco
1345               sis(il,ic)=sisca
1346               asy(il,ic)=asymf
1347               bre(il,ic)=refr
1348               bim(il,ic)=refi
1349            end if   
1350         end do
1351         read (ntap,'(7(/))')
1352         it=1
1353         ende=.false.
1354         do while (.not.ende)
1355            read (ntap,510,end=511)
1356     *            thet,(pha(il,ic,it),ila=1,ilamb)
1357  510       format(e11.3,1x,70e10.3)
1358            it=it+1
1359         end do
1360  511    ntheta=it-1
1361 
1362c ENDE NEUER INPUT 
1363 
1364            bac(il,ic)=pha(il,ic,ntheta)
1365           
1366            extbuf(ibuf)=ext(il,ic)           
1367            scabuf(ibuf)=sca(il,ic)           
1368            absbuf(ibuf)=abs(il,ic)           
1369            sisbuf(ibuf)=sis(il,ic)           
1370            asybuf(ibuf)=asy(il,ic)           
1371            brebuf(ibuf)=bre(il,ic)           
1372            bimbuf(ibuf)=bim(il,ic)           
1373            bacbuf(ibuf)=bac(il,ic)           
1374                   
1375            close (ntap)
1376           
1377c            print*,'closed file ',ntap
1378                       
1379         end if
1380         end do
1381      end do
1382
1383ccccc -----------------------------------------------------------------c
1384c     Formate                                                          c
1385ccccc -----------------------------------------------------------------c
1386
1387  100  format(8e10.3)
1388  200  format(a1)
1389  220  format(a10)
1390  300  format(15x,f6.3,/,8x,e10.4,7x,e10.4,7x,e10.4,5x,f7.4,7x,f6.4,/)
1391  301  format(8x,f6.3//)
1392  303  format(7x,e10.3,7x,e10.3,7x,e10.3,7x,f7.4/)
1393  400  format(12(/))
1394 1010  format(70X,e10.3)
1395
1396       return
1397       end
1398
1399CCCCC *****************************************************************C
1400      subroutine optpar (ilamb,ihum)
1401C     *****************************************************************C
1402C                                                                      c
1403C     -----------------------------------------------------------------C
1404C     Berechnung und Ausdruck der gewnschten optischen Parameter      C
1405CCCCC -----------------------------------------------------------------C
1406
1407      integer prnr,acnr,rht
1408      real n,numden
1409      REAL EXTN(2),ABSN(2),SCAN(2),PF18N(2),supf(112),phafu(112,2)
1410      REAL EXTA(2),ABSA(2),SCAA(2),SSA(2),ASF(2),PF18A(2)
1411      real scar(2),absr(2),omer(2)
1412      character*3  atn,pat
1413      character*4  comnam
1414      character*8  opanam,optnam
1415      character*20 catyp
1416      character*30 typnam
1417      common /atyp/   natyp,mcomp,ncomp(10),numden,
1418     *                catyp,typnam,comnam(20)
1419      common /norm/   norm,mixnor
1420      common /layer/  nlay,mlay,nltyp(10),parlay(10,2),boundl(10),
1421     *                boundu(10)
1422      common /profi/ nil(10),hfta(10),hstra(10),
1423     *                h0(2,10),h1(2,10),hm(2,10)
1424      common /opar/   mopar,jnopar(13),nop,opanam(13),optnam(13)
1425      common /wavel/  mlamb,alamb(61),niw
1426      common /mipoi/  latx,lonx,nl,prnr,rht(2),n(2),
1427     *                njc(2),acnr(5,2),acmr(5,2),nh(2),atn(2),pat(2)
1428      common /oppoi/  ext(2,5),sca(2,5),abs(2,5),sis(2,5),asy(2,5),
1429     *                bac(2,5),pha(2,5,112),bre(2,5),bim(2,5)
1430      COMMON /FTASTR/ EXTFTA(61),EXTSTR(61)
1431      COMMON /TEST/   ITEST(2)
1432      common /out/    oparam(10,2),phaf(112,2)
1433      common /prog/   nprog
1434      common /angle/  jnangle(112),angle(112),nia
1435      common /masse/  smas(10,8),smag(8)
1436
1437CCCCC ------------------------------------------------------------------C
1438C     MISCHEN DES AEROSOL-TYPS                                          C
1439C     SUMM(E,A,S) : SUMME EXTINCTION, ABSORPTION, SCATTERING            C
1440C     SUPF18      : SUMME DES RUECKSTREUKOEFFIZIENTEN                   C
1441C     SUMASF      : ZWISCHENSUMME DES ASYMMETRIEFAKTORS (ASF)           C
1442C     SUMASF      : ZWISCHENSUMME DER SINGLE SCAT. ALB. (SSA)           C
1443CCCCC ------------------------------------------------------------------C
1444
1445
1446      DO 10 L=1,NL
1447
1448      SUMME = 0.
1449      SUMMA = 0.
1450      SUMMS = 0.
1451      SUMSSA = 0.
1452      SUMASF = 0.
1453      SUPF18 = 0.
1454      if (jnopar(10).eq.1) then
1455      do it=1,112
1456         supf(it)=0.
1457      end do
1458      end if
1459
1460      DO 20 JC=1,NJC(L)
1461
1462c      print*,' Berechnung der Summen'
1463
1464      SUMME = SUMME + ACMR(JC,L)*EXT(l,jc)
1465      SUMMA = SUMMA + ACMR(JC,L)*ABS(l,jc)
1466      SUMMS = SUMMS + ACMR(JC,L)*SCA(l,jc)
1467      SUMSSA = SUMSSA + ACMR(JC,L)*sis(l,jc)
1468     +                  *EXT(l,jc)
1469      SUMASF = SUMASF + ACMR(JC,L)*asy(l,jc)
1470     +                  *SCA(l,jc)
1471      SUPF18 = SUPF18 + ACMR(JC,L)*bac(l,jc)
1472      if (jnopar(10).eq.1) then
1473         do it=1,112
1474            supf(it)=supf(it)+acmr(jc,l)*pha(l,jc,it)
1475         end do
1476      end if
1477     
1478c      print*,jc,l,njc(l),summe,acmr(jc,l),ext(l,jc)
1479
1480   20 CONTINUE
1481
1482CCCCC -----------------------------------------------------------------C
1483C     Normierte  optische Parameter                                    c
1484CCCCC -----------------------------------------------------------------C
1485
1486c      print*,' Berechnung der normierten Werte'
1487
1488      EXTN(L) = SUMME
1489      ABSN(L) = SUMMA
1490      SCAN(L) = SUMMS
1491      PF18N(L) = SUPF18
1492      if (jnopar(10).eq.1) then
1493      do it=1,112
1494         phafu(it,l)=supf(it)
1495      end do
1496      end if
1497
1498      SSA(L) = SUMSSA/SUMME
1499      ASF(L) = SUMASF/SUMMS
1500     
1501CCCCC -----------------------------------------------------------------C
1502C     ABSOLUTE OPTISCHE PARAMETER                                      C
1503CCCCC -----------------------------------------------------------------C
1504
1505c      print*,' Berechnung der absoluten Werte'
1506
1507      EXTA(L)= EXTN(L) * N(L)
1508      ABSA(L)= ABSN(L) * N(L)
1509      SCAA(L)= SCAN(L) * N(L)
1510      PF18A(L) = PF18N(L)* N(L)
1511      if (jnopar(10).eq.1.and.norm.eq.1) then
1512      do it=1,112
1513         phafu(it,l)=phafu(it,l)*n(l)
1514      end do
1515      end if
1516      if (norm.eq.1) then
1517         EXTN(L)= EXTA(L)
1518         ABSN(L)= ABSA(L)
1519         SCAN(L)= SCAA(L)
1520         PF18N(L) = PF18A(L)
1521      end if
1522
1523      if (jnopar(10).eq.1) then
1524         itp=1
1525         do it=1,112
1526            if (jnangle(it).eq.1) then
1527               phaf(itp,l)=phafu(it,l)
1528               itp=itp+1
1529            end if
1530         end do
1531      end if
1532
1533      if (jnopar(11).eq.1) then
1534         scar(l)=scaa(l)/smag(ihum)*1000.  ! Einheit m**2/g
1535      end if
1536
1537      if (jnopar(12).eq.1) then
1538         absr(l)=absa(l)/smag(ihum)*1000.
1539      end if
1540
1541      if (jnopar(13).eq.1) then
1542         kc=0
1543         do jc=1,njc(l)
1544            if (ncomp(jc).eq.3) kc=jc
1545         end do
1546         if (kc.ne.0) then
1547            omer(l)=ssa(l)/smas(kc,ihum)
1548         else
1549            omer(l)=99.
1550         end if
1551      end if
1552
1553CCCCC -----------------------------------------------------------------C
1554C     AUSGABE DER DATEN                                                C
1555CCCCC -----------------------------------------------------------------C
1556
1557      iop=0
1558      kop=0
1559      if (jnopar(1).eq.1) then
1560         iop=iop+1
1561         oparam(iop,l)=extn(l)
1562      end if
1563      if (jnopar(2).eq.1) then
1564         iop=iop+1
1565         oparam(iop,l)=scan(l)
1566      end if
1567      if (jnopar(3).eq.1) then
1568         iop=iop+1
1569         oparam(iop,l)=absn(l)
1570      end if
1571      if (jnopar(4).eq.1) then
1572         iop=iop+1
1573         oparam(iop,l)=ssa(l)
1574      end if
1575      if (jnopar(5).eq.1) then
1576         iop=iop+1
1577         oparam(iop,l)=asf(l)
1578      end if
1579      if (jnopar(9).eq.1) then
1580         iop=iop+1
1581         if (jnopar(6).eq.1) kop=kop+1
1582         if (jnopar(7).eq.1) kop=kop+1
1583         if (jnopar(8).eq.1) kop=kop+1
1584         kop=kop+iop
1585         oparam(kop,l)=exta(l)/pf18a(l)
1586      end if
1587      if (jnopar(11).eq.1) then
1588         iop=iop+1
1589         oparam(iop,l)=scar(l)
1590      end if
1591      if (jnopar(12).eq.1) then
1592         iop=iop+1
1593         oparam(iop,l)=absr(l)
1594      end if
1595      if (jnopar(13).eq.1) then
1596         iop=iop+1
1597         oparam(iop,l)=omer(l)
1598      end if
1599
1600   10 CONTINUE
1601
1602CCCCC -----------------------------------------------------------------C
1603C     OPTISCHE DICKE                                                   C
1604CCCCC -----------------------------------------------------------------C
1605
1606      if (jnopar(6).eq.1) then
1607
1608CCCCC -----------------------------------------------------------------C
1609c     Bestimmung von HM, HFTA, HSTR, EXTFTA, EXTSTR aus den            c
1610c     eingelesenen Werten in /layer/ fuer RAWOPT                       c
1611CCCCC -----------------------------------------------------------------C
1612
1613      if (nprog.eq.2) then
1614         do il=1,nl
1615            if (nltyp(il).eq.1) then
1616               hm(il,1)=parlay(il,1)
1617            else if (nltyp(il).eq.2) then
1618               hm(il,1) = parlay(il,2)*
1619     *                    (exp(-boundl(il)/parlay(il,2))+
1620     *                     exp(-boundu(il)/parlay(il,2)))
1621            end if
1622         end do
1623         hfta(1)=boundu(nlay-2)-boundl(nlay-2)
1624         hstra(1)=boundu(nlay-1)-boundl(nlay-1)
1625         extfta(ilamb)=parlay(nlay-2,1)
1626         extstr(ilamb)=parlay(nlay-1,1)
1627      end if
1628     
1629      hu = h1(nl,prnr)
1630      ho = hu + hfta(prnr)
1631      z = 8.
1632      hftae = z * ( exp(-hu/z) - exp(-ho/z) )
1633
1634      ODEPTH = 0.
1635
1636      DO IL=1,nl
1637         ODEPTH = ODEPTH + EXTA(IL) * HM(IL,prnr)
1638      end do
1639
1640CCCCC -----------------------------------------------------------------C
1641C                     + FREE TROP. AEROSOL                             C
1642CCCCC -----------------------------------------------------------------C
1643
1644         ODEPTH = ODEPTH + EXTFTA(ilamb)*HFTAE     
1645     +                   + EXTSTR(ilamb)*HSTRA(prnr)
1646
1647c        do il=1,nl
1648c        print*,'   exta= ',exta(il),'    hm(il)= ',hm(il,prnr)
1649c        end do
1650c        print*,' ilamb= ' ,ilamb
1651c        print*,' extfta= ',extfta(ilamb),' hftae= ',hftae
1652c        print*,' extstr= ',extstr(ilamb),' hstr= ',hstra(prnr)
1653
1654      odeptha=odepth/alog(10.)
1655
1656      turbr=0.008569*alamb(ilamb)**(-4)*(1.+0.0113*alamb(ilamb)**
1657     *         (-2)+0.00013*alamb(ilamb)**(-4))
1658
1659      turbf=(odepth+turbr)/turbr
1660
1661      if (jnopar(9).eq.1) then
1662         kop=iop
1663      else
1664         kop=iop+1
1665c        print*,' exta= ', exta(il),' hm= ',hm(il,1)
1666      end if
1667
1668      if (jnopar(6).eq.1) then
1669         oparam(kop,1)=odepth
1670         oparam(kop,2)=0.
1671         kop=kop+1
1672         iop=iop+1
1673      end if
1674
1675      if (jnopar(7).eq.1) then
1676         oparam(kop,1)=odeptha
1677         oparam(kop,2)=0.
1678         kop=kop+1
1679         iop=iop+1
1680      end if
1681
1682      if (jnopar(8).eq.1) then
1683         oparam(kop,1)=turbf
1684         oparam(kop,2)=0.
1685         iop=iop+1
1686      end if
1687      end if
1688
1689c      print*,'Aufruf von out4'
1690
1691      call out4(iop,ihum)
1692
1693c      print*,'Ende out4' 
1694
1695      RETURN
1696      END
1697
1698ccccc *****************************************************************c
1699      subroutine out4(iop,ihum)
1700c     *****************************************************************c
1701c                                                                            c
1702c     -----------------------------------------------------------------c
1703c     AUSGABE DER DATEN fr atlopt                                               c
1704ccccc -----------------------------------------------------------------c
1705
1706      integer prnr,acnr,njc,rht
1707      real n
1708
1709      character*2 chum
1710      character opanam*8,atn*3,pat*3,optnam*8
1711
1712      common /angle/  jnangle(112),angle(112),nia
1713      common /wavel/  mlamb,alamb(61),niw
1714      common /hum/    khum(8),ahum(8),nih,nhum(8),mhum,chum(8)
1715      common /out/    oparam(10,2),phaf(112,2)
1716      common /opar/   mopar,jnopar(13),nop,opanam(13),optnam(13)
1717      common /mipoi/  latx,lonx,nl,prnr,rht(2),n(2),
1718     *                njc(2),acnr(5,2),acmr(5,2),nh(2),atn(2),pat(2)
1719
1720
1721      do l=1,nl
1722         if (nih.ne.0) rht(l)=ahum(ihum)
1723         if (nop.gt.iop) then
1724            if (jnopar(9).eq.1.and.jnopar(10).eq.1) then
1725               oparam(iop,l)=oparam(nop-1,l)
1726            else if (jnopar(9).eq.1.and.jnopar(10).eq.0) then
1727               oparam(iop,l)=oparam(nop,l)
1728            end if
1729         end if
1730
1731         if (iop.le.10) then
1732            if (l.eq.1) then
1733               write (10,2020) latx,lonx,nl,
1734     *                         (oparam(ip,l),ip=1,iop)
1735 2020          FORMAT(2(1x,I4),i3,1p3e10.3,0p3e10.3,1pe10.3)
1736            else
1737               write (10,3020) 
1738     *                         (oparam(ip,l),ip=1,iop)
1739 3020          FORMAT(13x,1p3e10.3,0p3e10.3,1pe10.3)
1740            end if
1741         else
1742            if (l.eq.1) then
1743               write (10,2020) latx,lonx,nl,
1744     *                         (oparam(ip,l),ip=1,5)
1745               write (10,2030) (oparam(ip,l),ip=6,iop)
1746 2030          FORMAT(11x,1p10e10.3)
1747            else
1748               write (10,3040) 
1749     *                         (oparam(ip,l),ip=1,5)
1750 3040          FORMAT(13x,10e10.3)
1751               write (10,2030) (oparam(ip,l),ip=6,iop)
1752            end if
1753         end if
1754
1755         if(jnopar(10).eq.1) then
1756            write(10,4002)
1757 4002       format('   phase function')
1758            write(10,4010) (phaf(it,l),it=1,nia)
1759 4010       format(8e10.3)
1760            write(10,*) ' '
1761         end if
1762      end do
1763
1764      return
1765      end
Note: See TracBrowser for help on using the repository browser.