source: LMDZ5/trunk/libf/phylmd/cosp/icarus.F @ 2654

Last change on this file since 2654 was 2428, checked in by idelkadi, 9 years ago

Mise a jour du simulateur COSP (passage de la version v3.2 a la version v1.4) :

  • mise a jour des sources pour ISCCP, CALIPSO et PARASOL
  • prise en compte des changements de phases pour les nuages (Calipso)
  • rajout de plusieurs diagnostiques (fraction nuageuse en fonction de la temperature, ...)

http://lmdz.lmd.jussieu.fr/Members/aidelkadi/cosp

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 42.9 KB
Line 
1! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $
2! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/icarus-scops-4.1-bsd/icarus.f $
3      SUBROUTINE ICARUS(
4     &     debug,
5     &     debugcol,
6     &     npoints,
7     &     sunlit,
8     &     nlev,
9     &     ncol,
10     &     pfull,
11     &     phalf,
12     &     qv,
13     &     cc,
14     &     conv,
15     &     dtau_s,
16     &     dtau_c,
17     &     top_height,
18     &     top_height_direction,
19     &     overlap,
20     &     frac_out,
21     &     skt,
22     &     emsfc_lw,
23     &     at,
24     &     dem_s,
25     &     dem_c,
26     &     fq_isccp,
27     &     totalcldarea,
28     &     meanptop,
29     &     meantaucld,
30     &     meanalbedocld,
31     &     meantb,
32     &     meantbclr,
33     &     boxtau,
34     &     boxptop
35     &)
36
37!$Id: icarus.f,v 4.1 2010/05/27 16:30:18 hadmw Exp $
38
39! *****************************COPYRIGHT****************************
40! (c) 2009, Lawrence Livermore National Security Limited Liability
41! Corporation.
42! All rights reserved.
43!
44! Redistribution and use in source and binary forms, with or without
45! modification, are permitted provided that the
46! following conditions are met:
47!
48!     * Redistributions of source code must retain the above
49!       copyright  notice, this list of conditions and the following
50!       disclaimer.
51!     * Redistributions in binary form must reproduce the above
52!       copyright notice, this list of conditions and the following
53!       disclaimer in the documentation and/or other materials
54!       provided with the distribution.
55!     * Neither the name of the Lawrence Livermore National Security
56!       Limited Liability Corporation nor the names of its
57!       contributors may be used to endorse or promote products
58!       derived from this software without specific prior written
59!       permission.
60!
61! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
62! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
63! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
64! A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
65! OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
66! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
67! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
68! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
69! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
70! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
71! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 
72!
73! *****************************COPYRIGHT*******************************
74! *****************************COPYRIGHT*******************************
75! *****************************COPYRIGHT*******************************
76! *****************************COPYRIGHT*******************************
77
78      implicit none
79
80!     NOTE:   the maximum number of levels and columns is set by
81!             the following parameter statement
82
83      INTEGER ncolprint
84     
85!     -----
86!     Input
87!     -----
88
89      INTEGER npoints       !  number of model points in the horizontal
90      INTEGER nlev          !  number of model levels in column
91      INTEGER ncol          !  number of subcolumns
92
93      INTEGER sunlit(npoints) !  1 for day points, 0 for night time
94
95      REAL pfull(npoints,nlev)
96                       !  pressure of full model levels (Pascals)
97                  !  pfull(npoints,1) is top level of model
98                  !  pfull(npoints,nlev) is bot of model
99
100      REAL phalf(npoints,nlev+1)
101                  !  pressure of half model levels (Pascals)
102                  !  phalf(npoints,1) is top of model
103                  !  phalf(npoints,nlev+1) is the surface pressure
104
105      REAL qv(npoints,nlev)
106                  !  water vapor specific humidity (kg vapor/ kg air)
107                  !         on full model levels
108
109      REAL cc(npoints,nlev)   
110                  !  input cloud cover in each model level (fraction)
111                  !  NOTE:  This is the HORIZONTAL area of each
112                  !         grid box covered by clouds
113
114      REAL conv(npoints,nlev)
115                  !  input convective cloud cover in each model
116                  !   level (fraction)
117                  !  NOTE:  This is the HORIZONTAL area of each
118                  !         grid box covered by convective clouds
119
120      REAL dtau_s(npoints,nlev)
121                  !  mean 0.67 micron optical depth of stratiform
122                !  clouds in each model level
123                  !  NOTE:  this the cloud optical depth of only the
124                  !  cloudy part of the grid box, it is not weighted
125                  !  with the 0 cloud optical depth of the clear
126                  !         part of the grid box
127
128      REAL dtau_c(npoints,nlev)
129                  !  mean 0.67 micron optical depth of convective
130                !  clouds in each
131                  !  model level.  Same note applies as in dtau_s.
132
133      INTEGER overlap                   !  overlap type
134                              !  1=max
135                              !  2=rand
136                              !  3=max/rand
137
138      INTEGER top_height                !  1 = adjust top height using both a computed
139                                        !  infrared brightness temperature and the visible
140                              !  optical depth to adjust cloud top pressure. Note
141                              !  that this calculation is most appropriate to compare
142                              !  to ISCCP data during sunlit hours.
143                                        !  2 = do not adjust top height, that is cloud top
144                                        !  pressure is the actual cloud top pressure
145                                        !  in the model
146                              !  3 = adjust top height using only the computed
147                              !  infrared brightness temperature. Note that this
148                              !  calculation is most appropriate to compare to ISCCP
149                              !  IR only algortihm (i.e. you can compare to nighttime
150                              !  ISCCP data with this option)
151
152      INTEGER top_height_direction ! direction for finding atmosphere pressure level
153                                 ! with interpolated temperature equal to the radiance
154                                 ! determined cloud-top temperature
155                                 !
156                                 ! 1 = find the *lowest* altitude (highest pressure) level
157                                 ! with interpolated temperature equal to the radiance
158                                 ! determined cloud-top temperature
159                                 !
160                                 ! 2 = find the *highest* altitude (lowest pressure) level
161                                 ! with interpolated temperature equal to the radiance
162                                 ! determined cloud-top temperature
163                                 !
164                                 ! ONLY APPLICABLE IF top_height EQUALS 1 or 3
165                                 !                               !
166                                 ! 1 = old setting: matches all versions of
167                                 ! ISCCP simulator with versions numbers 3.5.1 and lower
168                                 !
169                                 ! 2 = default setting: for version numbers 4.0 and higher
170!
171!     The following input variables are used only if top_height = 1 or top_height = 3
172!
173      REAL skt(npoints)                 !  skin Temperature (K)
174      REAL emsfc_lw                     !  10.5 micron emissivity of surface (fraction)                                           
175      REAL at(npoints,nlev)                   !  temperature in each model level (K)
176      REAL dem_s(npoints,nlev)                !  10.5 micron longwave emissivity of stratiform
177                              !  clouds in each
178                                        !  model level.  Same note applies as in dtau_s.
179      REAL dem_c(npoints,nlev)                  !  10.5 micron longwave emissivity of convective
180                              !  clouds in each
181                                        !  model level.  Same note applies as in dtau_s.
182
183      REAL frac_out(npoints,ncol,nlev) ! boxes gridbox divided up into
184                              ! Equivalent of BOX in original version, but
185                              ! indexed by column then row, rather than
186                              ! by row then column
187
188
189
190!     ------
191!     Output
192!     ------
193
194      REAL fq_isccp(npoints,7,7)        !  the fraction of the model grid box covered by
195                                        !  each of the 49 ISCCP D level cloud types
196
197      REAL totalcldarea(npoints)        !  the fraction of model grid box columns
198                                        !  with cloud somewhere in them.  NOTE: This diagnostic
199                                        ! does not count model clouds with tau < isccp_taumin
200                              ! Thus this diagnostic does not equal the sum over all entries of fq_isccp.
201                              ! However, this diagnostic does equal the sum over entries of fq_isccp with
202                              ! itau = 2:7 (omitting itau = 1)
203     
204     
205      ! The following three means are averages only over the cloudy areas with tau > isccp_taumin. 
206      ! If no clouds with tau > isccp_taumin are in grid box all three quantities should equal zero.     
207                             
208      REAL meanptop(npoints)            !  mean cloud top pressure (mb) - linear averaging
209                                        !  in cloud top pressure.
210                             
211      REAL meantaucld(npoints)          !  mean optical thickness
212                                        !  linear averaging in albedo performed.
213     
214      real meanalbedocld(npoints)        ! mean cloud albedo
215                                        ! linear averaging in albedo performed
216                                       
217      real meantb(npoints)              ! mean all-sky 10.5 micron brightness temperature
218     
219      real meantbclr(npoints)           ! mean clear-sky 10.5 micron brightness temperature
220     
221      REAL boxtau(npoints,ncol)         !  optical thickness in each column
222     
223      REAL boxptop(npoints,ncol)        !  cloud top pressure (mb) in each column
224                             
225                                                                                         
226!
227!     ------
228!     Working variables added when program updated to mimic Mark Webb's PV-Wave code
229!     ------
230
231      REAL dem(npoints,ncol),bb(npoints)     !  working variables for 10.5 micron longwave
232                              !  emissivity in part of
233                              !  gridbox under consideration
234
235      REAL ptrop(npoints)
236      REAL attrop(npoints)
237      REAL attropmin (npoints)
238      REAL atmax(npoints)
239      REAL btcmin(npoints)
240      REAL transmax(npoints)
241
242      INTEGER i,j,ilev,ibox,itrop(npoints)
243      INTEGER ipres(npoints)
244      INTEGER itau(npoints),ilev2
245      INTEGER acc(nlev,ncol)
246      INTEGER match(npoints,nlev-1)
247      INTEGER nmatch(npoints)
248      INTEGER levmatch(npoints,ncol)
249     
250      !variables needed for water vapor continuum absorption
251      real fluxtop_clrsky(npoints),trans_layers_above_clrsky(npoints)
252      real taumin(npoints)
253      real dem_wv(npoints,nlev), wtmair, wtmh20, Navo, grav, pstd, t0
254      real press(npoints), dpress(npoints), atmden(npoints)
255      real rvh20(npoints), wk(npoints), rhoave(npoints)
256      real rh20s(npoints), rfrgn(npoints)
257      real tmpexp(npoints),tauwv(npoints)
258     
259      character*1 cchar(6),cchar_realtops(6)
260      integer icycle
261      REAL tau(npoints,ncol)
262      LOGICAL box_cloudy(npoints,ncol)
263      REAL tb(npoints,ncol)
264      REAL ptop(npoints,ncol)
265      REAL emcld(npoints,ncol)
266      REAL fluxtop(npoints,ncol)
267      REAL trans_layers_above(npoints,ncol)
268      real isccp_taumin,fluxtopinit(npoints),tauir(npoints)
269      REAL albedocld(npoints,ncol)
270      real boxarea
271      integer debug       ! set to non-zero value to print out inputs
272                    ! with step debug
273      integer debugcol    ! set to non-zero value to print out column
274                    ! decomposition with step debugcol
275      integer rangevec(npoints),rangeerror
276
277      integer index1(npoints),num1,jj,k1,k2
278      real rec2p13,tauchk,logp,logp1,logp2,atd
279      real output_missing_value
280
281      character*10 ftn09
282     
283      DATA isccp_taumin / 0.3 /
284      DATA output_missing_value / -1.E+30 /
285      DATA cchar / ' ','-','1','+','I','+'/
286      DATA cchar_realtops / ' ',' ','1','1','I','I'/
287
288!     ------ End duplicate definitions common to wrapper routine
289
290      tauchk = -1.*log(0.9999999)
291      rec2p13=1./2.13
292
293      ncolprint=0
294
295      if ( debug.ne.0 ) then
296          j=1
297          write(6,'(a10)') 'j='
298          write(6,'(8I10)') j
299          write(6,'(a10)') 'debug='
300          write(6,'(8I10)') debug
301          write(6,'(a10)') 'debugcol='
302          write(6,'(8I10)') debugcol
303          write(6,'(a10)') 'npoints='
304          write(6,'(8I10)') npoints
305          write(6,'(a10)') 'nlev='
306          write(6,'(8I10)') nlev
307          write(6,'(a10)') 'ncol='
308          write(6,'(8I10)') ncol
309          write(6,'(a11)') 'top_height='
310          write(6,'(8I10)') top_height
311          write(6,'(a21)') 'top_height_direction='
312          write(6,'(8I10)') top_height_direction
313          write(6,'(a10)') 'overlap='
314          write(6,'(8I10)') overlap
315          write(6,'(a10)') 'emsfc_lw='
316          write(6,'(8f10.2)') emsfc_lw
317        do j=1,npoints,debug
318          write(6,'(a10)') 'j='
319          write(6,'(8I10)') j
320          write(6,'(a10)') 'sunlit='
321          write(6,'(8I10)') sunlit(j)
322          write(6,'(a10)') 'pfull='
323          write(6,'(8f10.2)') (pfull(j,i),i=1,nlev)
324          write(6,'(a10)') 'phalf='
325          write(6,'(8f10.2)') (phalf(j,i),i=1,nlev+1)
326          write(6,'(a10)') 'qv='
327          write(6,'(8f10.3)') (qv(j,i),i=1,nlev)
328          write(6,'(a10)') 'cc='
329          write(6,'(8f10.3)') (cc(j,i),i=1,nlev)
330          write(6,'(a10)') 'conv='
331          write(6,'(8f10.2)') (conv(j,i),i=1,nlev)
332          write(6,'(a10)') 'dtau_s='
333          write(6,'(8g12.5)') (dtau_s(j,i),i=1,nlev)
334          write(6,'(a10)') 'dtau_c='
335          write(6,'(8f10.2)') (dtau_c(j,i),i=1,nlev)
336          write(6,'(a10)') 'skt='
337          write(6,'(8f10.2)') skt(j)
338          write(6,'(a10)') 'at='
339          write(6,'(8f10.2)') (at(j,i),i=1,nlev)
340          write(6,'(a10)') 'dem_s='
341          write(6,'(8f10.3)') (dem_s(j,i),i=1,nlev)
342          write(6,'(a10)') 'dem_c='
343          write(6,'(8f10.3)') (dem_c(j,i),i=1,nlev)
344        enddo
345      endif
346
347!     ---------------------------------------------------!
348
349      if (ncolprint.ne.0) then
350      do j=1,npoints,1000
351        write(6,'(a10)') 'j='
352        write(6,'(8I10)') j
353      enddo
354      endif
355
356      if (top_height .eq. 1 .or. top_height .eq. 3) then
357
358      do j=1,npoints
359          ptrop(j)=5000.
360          attropmin(j) = 400.
361          atmax(j) = 0.
362          attrop(j) = 120.
363          itrop(j) = 1
364      enddo
365
366      do 12 ilev=1,nlev
367        do j=1,npoints
368         if (pfull(j,ilev) .lt. 40000. .and.
369     &          pfull(j,ilev) .gt.  5000. .and.
370     &          at(j,ilev) .lt. attropmin(j)) then
371                ptrop(j) = pfull(j,ilev)
372                attropmin(j) = at(j,ilev)
373                attrop(j) = attropmin(j)
374                itrop(j)=ilev
375           end if
376        enddo
37712    continue
378
379      do 13 ilev=1,nlev
380        do j=1,npoints
381           if (at(j,ilev) .gt. atmax(j) .and.
382     &              ilev  .ge. itrop(j)) atmax(j)=at(j,ilev)
383        enddo
38413    continue
385
386      end if
387
388
389      if (top_height .eq. 1 .or. top_height .eq. 3) then
390          do j=1,npoints
391              meantb(j) = 0.
392              meantbclr(j) = 0.
393          end do
394      else
395          do j=1,npoints
396              meantb(j) = output_missing_value
397              meantbclr(j) = output_missing_value
398          end do
399      end if
400     
401!     -----------------------------------------------------!
402
403!     ---------------------------------------------------!
404
405      do ilev=1,nlev
406        do j=1,npoints
407
408          rangevec(j)=0
409
410          if (cc(j,ilev) .lt. 0. .or. cc(j,ilev) .gt. 1.) then
411!           error = cloud fraction less than zero
412!           error = cloud fraction greater than 1
413            rangevec(j)=rangevec(j)+1
414          endif
415
416          if (conv(j,ilev) .lt. 0. .or. conv(j,ilev) .gt. 1.) then
417!           ' error = convective cloud fraction less than zero'
418!           ' error = convective cloud fraction greater than 1'
419            rangevec(j)=rangevec(j)+2
420          endif
421
422          if (dtau_s(j,ilev) .lt. 0.) then
423!           ' error = stratiform cloud opt. depth less than zero'
424            rangevec(j)=rangevec(j)+4
425          endif
426
427          if (dtau_c(j,ilev) .lt. 0.) then
428!           ' error = convective cloud opt. depth less than zero'
429            rangevec(j)=rangevec(j)+8
430          endif
431
432          if (dem_s(j,ilev) .lt. 0. .or. dem_s(j,ilev) .gt. 1.) then
433!             ' error = stratiform cloud emissivity less than zero'
434!             ' error = stratiform cloud emissivity greater than 1'
435            rangevec(j)=rangevec(j)+16
436          endif
437
438          if (dem_c(j,ilev) .lt. 0. .or. dem_c(j,ilev) .gt. 1.) then
439!             ' error = convective cloud emissivity less than zero'
440!             ' error = convective cloud emissivity greater than 1'
441              rangevec(j)=rangevec(j)+32
442          endif
443        enddo
444
445        rangeerror=0
446        do j=1,npoints
447            rangeerror=rangeerror+rangevec(j)
448        enddo
449
450        if (rangeerror.ne.0) then
451              write (6,*) 'Input variable out of range'
452              write (6,*) 'rangevec:'
453              write (6,*) rangevec
454              STOP
455        endif
456      enddo
457
458!
459!     ---------------------------------------------------!
460
461     
462!
463!     ---------------------------------------------------!
464!     COMPUTE CLOUD OPTICAL DEPTH FOR EACH COLUMN and
465!     put into vector tau
466 
467      !initialize tau and albedocld to zero
468      do 15 ibox=1,ncol
469        do j=1,npoints
470            tau(j,ibox)=0.
471          albedocld(j,ibox)=0.
472          boxtau(j,ibox)=output_missing_value
473          boxptop(j,ibox)=output_missing_value
474          box_cloudy(j,ibox)=.false.
475        enddo
47615    continue
477
478      !compute total cloud optical depth for each column     
479      do ilev=1,nlev
480            !increment tau for each of the boxes
481            do ibox=1,ncol
482              do j=1,npoints
483                 if (frac_out(j,ibox,ilev).eq.1) then
484                        tau(j,ibox)=tau(j,ibox)
485     &                     + dtau_s(j,ilev)
486                 endif
487                 if (frac_out(j,ibox,ilev).eq.2) then
488                        tau(j,ibox)=tau(j,ibox)
489     &                     + dtau_c(j,ilev)
490                 end if
491              enddo
492            enddo ! ibox
493      enddo ! ilev
494          if (ncolprint.ne.0) then
495
496              do j=1,npoints ,1000
497                write(6,'(a10)') 'j='
498                write(6,'(8I10)') j
499                write(6,'(i2,1X,8(f7.2,1X))')
500     &          ilev,
501     &          (tau(j,ibox),ibox=1,ncolprint)
502              enddo
503          endif
504!
505!     ---------------------------------------------------!
506
507
508
509!     
510!     ---------------------------------------------------!
511!     COMPUTE INFRARED BRIGHTNESS TEMPERUATRES
512!     AND CLOUD TOP TEMPERATURE SATELLITE SHOULD SEE
513!
514!     again this is only done if top_height = 1 or 3
515!
516!     fluxtop is the 10.5 micron radiance at the top of the
517!              atmosphere
518!     trans_layers_above is the total transmissivity in the layers
519!             above the current layer
520!     fluxtop_clrsky(j) and trans_layers_above_clrsky(j) are the clear
521!             sky versions of these quantities.
522
523      if (top_height .eq. 1 .or. top_height .eq. 3) then
524
525
526        !----------------------------------------------------------------------
527        !   
528        !             DO CLEAR SKY RADIANCE CALCULATION FIRST
529        !
530        !compute water vapor continuum emissivity
531        !this treatment follows Schwarkzopf and Ramasamy
532        !JGR 1999,vol 104, pages 9467-9499.
533        !the emissivity is calculated at a wavenumber of 955 cm-1,
534        !or 10.47 microns
535        wtmair = 28.9644
536        wtmh20 = 18.01534
537        Navo = 6.023E+23
538        grav = 9.806650E+02
539        pstd = 1.013250E+06
540        t0 = 296.
541        if (ncolprint .ne. 0)
542     &         write(6,*)  'ilev   pw (kg/m2)   tauwv(j)      dem_wv'
543        do 125 ilev=1,nlev
544          do j=1,npoints
545               !press and dpress are dyne/cm2 = Pascals *10
546               press(j) = pfull(j,ilev)*10.
547               dpress(j) = (phalf(j,ilev+1)-phalf(j,ilev))*10
548               !atmden = g/cm2 = kg/m2 / 10
549               atmden(j) = dpress(j)/grav
550               rvh20(j) = qv(j,ilev)*wtmair/wtmh20
551               wk(j) = rvh20(j)*Navo*atmden(j)/wtmair
552               rhoave(j) = (press(j)/pstd)*(t0/at(j,ilev))
553               rh20s(j) = rvh20(j)*rhoave(j)
554               rfrgn(j) = rhoave(j)-rh20s(j)
555               tmpexp(j) = exp(-0.02*(at(j,ilev)-t0))
556               tauwv(j) = wk(j)*1.e-20*(
557     &           (0.0224697*rh20s(j)*tmpexp(j)) +
558     &                (3.41817e-7*rfrgn(j)) )*0.98
559               dem_wv(j,ilev) = 1. - exp( -1. * tauwv(j))
560          enddo
561               if (ncolprint .ne. 0) then
562               do j=1,npoints ,1000
563               write(6,'(a10)') 'j='
564               write(6,'(8I10)') j
565               write(6,'(i2,1X,3(f8.3,3X))') ilev,
566     &           qv(j,ilev)*(phalf(j,ilev+1)-phalf(j,ilev))/(grav/100.),
567     &           tauwv(j),dem_wv(j,ilev)
568               enddo
569             endif
570125     continue
571
572        !initialize variables
573        do j=1,npoints
574          fluxtop_clrsky(j) = 0.
575          trans_layers_above_clrsky(j)=1.
576        enddo
577
578        do ilev=1,nlev
579          do j=1,npoints
580 
581            ! Black body emission at temperature of the layer
582
583              bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. )
584              !bb(j)= 5.67e-8*at(j,ilev)**4
585
586              ! increase TOA flux by flux emitted from layer
587              ! times total transmittance in layers above
588
589                fluxtop_clrsky(j) = fluxtop_clrsky(j)
590     &            + dem_wv(j,ilev)*bb(j)*trans_layers_above_clrsky(j)
591           
592                ! update trans_layers_above with transmissivity
593              ! from this layer for next time around loop
594
595                trans_layers_above_clrsky(j)=
596     &            trans_layers_above_clrsky(j)*(1.-dem_wv(j,ilev))
597                   
598
599          enddo   
600            if (ncolprint.ne.0) then
601             do j=1,npoints ,1000
602              write(6,'(a10)') 'j='
603              write(6,'(8I10)') j
604              write (6,'(a)') 'ilev:'
605              write (6,'(I2)') ilev
606   
607              write (6,'(a)')
608     &        'emiss_layer,100.*bb(j),100.*f,total_trans:'
609              write (6,'(4(f7.2,1X))') dem_wv(j,ilev),100.*bb(j),
610     &             100.*fluxtop_clrsky(j),trans_layers_above_clrsky(j)
611             enddo   
612            endif
613
614        enddo   !loop over level
615       
616        do j=1,npoints
617          !add in surface emission
618          bb(j)=1/( exp(1307.27/skt(j)) - 1. )
619          !bb(j)=5.67e-8*skt(j)**4
620
621          fluxtop_clrsky(j) = fluxtop_clrsky(j) + emsfc_lw * bb(j)
622     &     * trans_layers_above_clrsky(j)
623     
624          !clear sky brightness temperature
625          meantbclr(j) = 1307.27/(log(1.+(1./fluxtop_clrsky(j))))
626         
627        enddo
628
629        if (ncolprint.ne.0) then
630        do j=1,npoints ,1000
631          write(6,'(a10)') 'j='
632          write(6,'(8I10)') j
633          write (6,'(a)') 'id:'
634          write (6,'(a)') 'surface'
635
636          write (6,'(a)') 'emsfc,100.*bb(j),100.*f,total_trans:'
637          write (6,'(5(f7.2,1X))') emsfc_lw,100.*bb(j),
638     &      100.*fluxtop_clrsky(j),
639     &       trans_layers_above_clrsky(j), meantbclr(j)
640        enddo
641      endif
642   
643
644        !
645        !           END OF CLEAR SKY CALCULATION
646        !
647        !----------------------------------------------------------------
648
649
650
651        if (ncolprint.ne.0) then
652
653        do j=1,npoints ,1000
654            write(6,'(a10)') 'j='
655            write(6,'(8I10)') j
656            write (6,'(a)') 'ts:'
657            write (6,'(8f7.2)') (skt(j),ibox=1,ncolprint)
658   
659            write (6,'(a)') 'ta_rev:'
660            write (6,'(8f7.2)')
661     &       ((at(j,ilev2),ibox=1,ncolprint),ilev2=1,nlev)
662
663        enddo
664        endif
665        !loop over columns
666        do ibox=1,ncol
667          do j=1,npoints
668            fluxtop(j,ibox)=0.
669            trans_layers_above(j,ibox)=1.
670          enddo
671        enddo
672
673        do ilev=1,nlev
674              do j=1,npoints
675                ! Black body emission at temperature of the layer
676
677              bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. )
678              !bb(j)= 5.67e-8*at(j,ilev)**4
679              enddo
680
681            do ibox=1,ncol
682              do j=1,npoints
683
684              ! emissivity for point in this layer
685                if (frac_out(j,ibox,ilev).eq.1) then
686                dem(j,ibox)= 1. -
687     &             ( (1. - dem_wv(j,ilev)) * (1. -  dem_s(j,ilev)) )
688                else if (frac_out(j,ibox,ilev).eq.2) then
689                dem(j,ibox)= 1. -
690     &             ( (1. - dem_wv(j,ilev)) * (1. -  dem_c(j,ilev)) )
691                else
692                dem(j,ibox)=  dem_wv(j,ilev)
693                end if
694               
695
696                ! increase TOA flux by flux emitted from layer
697              ! times total transmittance in layers above
698
699                fluxtop(j,ibox) = fluxtop(j,ibox)
700     &            + dem(j,ibox) * bb(j)
701     &            * trans_layers_above(j,ibox)
702           
703                ! update trans_layers_above with transmissivity
704              ! from this layer for next time around loop
705
706                trans_layers_above(j,ibox)=
707     &            trans_layers_above(j,ibox)*(1.-dem(j,ibox))
708
709              enddo ! j
710            enddo ! ibox
711
712            if (ncolprint.ne.0) then
713              do j=1,npoints,1000
714              write (6,'(a)') 'ilev:'
715              write (6,'(I2)') ilev
716   
717              write(6,'(a10)') 'j='
718              write(6,'(8I10)') j
719              write (6,'(a)') 'emiss_layer:'
720              write (6,'(8f7.2)') (dem(j,ibox),ibox=1,ncolprint)
721       
722              write (6,'(a)') '100.*bb(j):'
723              write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint)
724       
725              write (6,'(a)') '100.*f:'
726              write (6,'(8f7.2)')
727     &         (100.*fluxtop(j,ibox),ibox=1,ncolprint)
728       
729              write (6,'(a)') 'total_trans:'
730              write (6,'(8f7.2)')
731     &          (trans_layers_above(j,ibox),ibox=1,ncolprint)
732            enddo
733          endif
734
735        enddo ! ilev
736
737
738          do j=1,npoints
739            !add in surface emission
740            bb(j)=1/( exp(1307.27/skt(j)) - 1. )
741            !bb(j)=5.67e-8*skt(j)**4
742          end do
743
744        do ibox=1,ncol
745          do j=1,npoints
746
747            !add in surface emission
748
749            fluxtop(j,ibox) = fluxtop(j,ibox)
750     &         + emsfc_lw * bb(j)
751     &         * trans_layers_above(j,ibox)
752           
753          end do
754        end do
755
756        !calculate mean infrared brightness temperature
757        do ibox=1,ncol
758          do j=1,npoints
759            meantb(j) = meantb(j)+1307.27/(log(1.+(1./fluxtop(j,ibox))))
760          end do
761        end do
762          do j=1, npoints
763            meantb(j) = meantb(j) / real(ncol)
764          end do       
765
766        if (ncolprint.ne.0) then
767
768          do j=1,npoints ,1000
769          write(6,'(a10)') 'j='
770          write(6,'(8I10)') j
771          write (6,'(a)') 'id:'
772          write (6,'(a)') 'surface'
773
774          write (6,'(a)') 'emiss_layer:'
775          write (6,'(8f7.2)') (dem(1,ibox),ibox=1,ncolprint)
776   
777          write (6,'(a)') '100.*bb(j):'
778          write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint)
779   
780          write (6,'(a)') '100.*f:'
781          write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint)
782         
783          write (6,'(a)') 'meantb(j):'
784          write (6,'(8f7.2)') (meantb(j),ibox=1,ncolprint)
785     
786          end do
787      endif
788   
789        !now that you have the top of atmosphere radiance account
790        !for ISCCP procedures to determine cloud top temperature
791
792        !account for partially transmitting cloud recompute flux
793        !ISCCP would see assuming a single layer cloud
794        !note choice here of 2.13, as it is primarily ice
795        !clouds which have partial emissivity and need the
796        !adjustment performed in this section
797        !
798      !If it turns out that the cloud brightness temperature
799      !is greater than 260K, then the liquid cloud conversion
800        !factor of 2.56 is used.
801      !
802        !Note that this is discussed on pages 85-87 of
803        !the ISCCP D level documentation (Rossow et al. 1996)
804           
805          do j=1,npoints 
806            !compute minimum brightness temperature and optical depth
807            btcmin(j) = 1. /  ( exp(1307.27/(attrop(j)-5.)) - 1. )
808          enddo
809        do ibox=1,ncol
810          do j=1,npoints 
811            transmax(j) = (fluxtop(j,ibox)-btcmin(j))
812     &                /(fluxtop_clrsky(j)-btcmin(j))
813          !note that the initial setting of tauir(j) is needed so that
814          !tauir(j) has a realistic value should the next if block be
815          !bypassed
816            tauir(j) = tau(j,ibox) * rec2p13
817            taumin(j) = -1. * log(max(min(transmax(j),0.9999999),0.001))
818
819          enddo
820
821          if (top_height .eq. 1) then
822            do j=1,npoints 
823              if (transmax(j) .gt. 0.001 .and.
824     &          transmax(j) .le. 0.9999999) then
825                fluxtopinit(j) = fluxtop(j,ibox)
826              tauir(j) = tau(j,ibox) *rec2p13
827              endif
828            enddo
829            do icycle=1,2
830              do j=1,npoints 
831                if (tau(j,ibox) .gt. (tauchk            )) then
832                if (transmax(j) .gt. 0.001 .and.
833     &            transmax(j) .le. 0.9999999) then
834                  emcld(j,ibox) = 1. - exp(-1. * tauir(j)  )
835                  fluxtop(j,ibox) = fluxtopinit(j) -   
836     &              ((1.-emcld(j,ibox))*fluxtop_clrsky(j))
837                  fluxtop(j,ibox)=max(1.E-06,
838     &              (fluxtop(j,ibox)/emcld(j,ibox)))
839                  tb(j,ibox)= 1307.27
840     &              / (log(1. + (1./fluxtop(j,ibox))))
841                  if (tb(j,ibox) .gt. 260.) then
842                  tauir(j) = tau(j,ibox) / 2.56
843                  end if                   
844                end if
845                end if
846              enddo
847            enddo
848               
849          endif
850       
851          do j=1,npoints
852            if (tau(j,ibox) .gt. (tauchk            )) then
853                !cloudy box
854                !NOTE: tb is the cloud-top temperature not infrared brightness temperature
855                !at this point in the code
856                tb(j,ibox)= 1307.27/ (log(1. + (1./fluxtop(j,ibox))))
857                if (top_height.eq.1.and.tauir(j).lt.taumin(j)) then
858                         tb(j,ibox) = attrop(j) - 5.
859                   tau(j,ibox) = 2.13*taumin(j)
860                end if
861            else
862                !clear sky brightness temperature
863                tb(j,ibox) = meantbclr(j)
864            end if
865          enddo ! j
866        enddo ! ibox
867
868        if (ncolprint.ne.0) then
869
870          do j=1,npoints,1000
871          write(6,'(a10)') 'j='
872          write(6,'(8I10)') j
873
874          write (6,'(a)') 'attrop:'
875          write (6,'(8f7.2)') (attrop(j))
876   
877          write (6,'(a)') 'btcmin:'
878          write (6,'(8f7.2)') (btcmin(j))
879   
880          write (6,'(a)') 'fluxtop_clrsky*100:'
881          write (6,'(8f7.2)')
882     &      (100.*fluxtop_clrsky(j))
883
884          write (6,'(a)') '100.*f_adj:'
885          write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint)
886   
887          write (6,'(a)') 'transmax:'
888          write (6,'(8f7.2)') (transmax(ibox),ibox=1,ncolprint)
889   
890          write (6,'(a)') 'tau:'
891          write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint)
892   
893          write (6,'(a)') 'emcld:'
894          write (6,'(8f7.2)') (emcld(j,ibox),ibox=1,ncolprint)
895   
896          write (6,'(a)') 'total_trans:'
897          write (6,'(8f7.2)')
898     &        (trans_layers_above(j,ibox),ibox=1,ncolprint)
899   
900          write (6,'(a)') 'total_emiss:'
901          write (6,'(8f7.2)')
902     &        (1.0-trans_layers_above(j,ibox),ibox=1,ncolprint)
903   
904          write (6,'(a)') 'total_trans:'
905          write (6,'(8f7.2)')
906     &        (trans_layers_above(j,ibox),ibox=1,ncolprint)
907   
908          write (6,'(a)') 'ppout:'
909          write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint)
910          enddo ! j
911      endif
912
913      end if
914
915!     ---------------------------------------------------!
916
917!     
918!     ---------------------------------------------------!
919!     DETERMINE CLOUD TOP PRESSURE
920!
921!     again the 2 methods differ according to whether
922!     or not you use the physical cloud top pressure (top_height = 2)
923!     or the radiatively determined cloud top pressure (top_height = 1 or 3)
924!
925
926      !compute cloud top pressure
927      do 30 ibox=1,ncol
928        !segregate according to optical thickness
929        if (top_height .eq. 1 .or. top_height .eq. 3) then 
930          !find level whose temperature
931          !most closely matches brightness temperature
932          do j=1,npoints
933            nmatch(j)=0
934          enddo
935          do 29 k1=1,nlev-1
936            if (top_height_direction .eq. 2) then
937              ilev = nlev - k1
938            else
939              ilev = k1
940            end if
941            !cdir nodep
942            do j=1,npoints
943             if (ilev .ge. itrop(j)) then
944              if ((at(j,ilev)   .ge. tb(j,ibox) .and.
945     &          at(j,ilev+1) .le. tb(j,ibox)) .or.
946     &          (at(j,ilev) .le. tb(j,ibox) .and.
947     &          at(j,ilev+1) .ge. tb(j,ibox))) then
948                nmatch(j)=nmatch(j)+1
949                match(j,nmatch(j))=ilev
950              end if 
951             end if                         
952            enddo
95329        continue
954
955          do j=1,npoints
956            if (nmatch(j) .ge. 1) then
957              k1 = match(j,nmatch(j))
958              k2 = k1 + 1
959              logp1 = log(pfull(j,k1))
960              logp2 = log(pfull(j,k2))
961              atd = max(tauchk,abs(at(j,k2) - at(j,k1)))
962              logp=logp1+(logp2-logp1)*abs(tb(j,ibox)-at(j,k1))/atd
963              ptop(j,ibox) = exp(logp)
964              if(abs(pfull(j,k1)-ptop(j,ibox)) .lt.
965     &            abs(pfull(j,k2)-ptop(j,ibox))) then
966                 levmatch(j,ibox)=k1
967              else
968                 levmatch(j,ibox)=k2
969              end if   
970            else
971              if (tb(j,ibox) .le. attrop(j)) then
972                ptop(j,ibox)=ptrop(j)
973                levmatch(j,ibox)=itrop(j)
974              end if
975              if (tb(j,ibox) .ge. atmax(j)) then
976                ptop(j,ibox)=pfull(j,nlev)
977                levmatch(j,ibox)=nlev
978              end if                               
979            end if
980          enddo ! j
981
982        else ! if (top_height .eq. 1 .or. top_height .eq. 3)
983 
984          do j=1,npoints     
985            ptop(j,ibox)=0.
986          enddo
987          do ilev=1,nlev
988            do j=1,npoints     
989              if ((ptop(j,ibox) .eq. 0. )
990     &           .and.(frac_out(j,ibox,ilev) .ne. 0)) then
991                ptop(j,ibox)=phalf(j,ilev)
992              levmatch(j,ibox)=ilev
993              end if
994            end do
995          end do
996        end if                           
997         
998        do j=1,npoints
999          if (tau(j,ibox) .le. (tauchk            )) then
1000            ptop(j,ibox)=0.
1001            levmatch(j,ibox)=0     
1002          endif
1003        enddo
1004
100530    continue
1006             
1007!
1008!
1009!     ---------------------------------------------------!
1010
1011
1012!     
1013!     ---------------------------------------------------!
1014!     DETERMINE ISCCP CLOUD TYPE FREQUENCIES
1015!
1016!     Now that ptop and tau have been determined,
1017!     determine amount of each of the 49 ISCCP cloud
1018!     types
1019!
1020!     Also compute grid box mean cloud top pressure and
1021!     optical thickness.  The mean cloud top pressure and
1022!     optical thickness are averages over the cloudy
1023!     area only. The mean cloud top pressure is a linear
1024!     average of the cloud top pressures.  The mean cloud
1025!     optical thickness is computed by converting optical
1026!     thickness to an albedo, averaging in albedo units,
1027!     then converting the average albedo back to a mean
1028!     optical thickness. 
1029!
1030
1031      !compute isccp frequencies
1032
1033      !reset frequencies
1034      do 38 ilev=1,7
1035      do 38 ilev2=1,7
1036        do j=1,npoints !
1037             if (sunlit(j).eq.1 .or. top_height .eq. 3) then
1038                fq_isccp(j,ilev,ilev2)= 0.
1039             else
1040                fq_isccp(j,ilev,ilev2)= output_missing_value
1041             end if
1042        enddo
104338    continue
1044
1045      !reset variables need for averaging cloud properties
1046      do j=1,npoints
1047        if (sunlit(j).eq.1 .or. top_height .eq. 3) then
1048             totalcldarea(j) = 0.
1049             meanalbedocld(j) = 0.
1050             meanptop(j) = 0.
1051             meantaucld(j) = 0.
1052        else
1053             totalcldarea(j) = output_missing_value
1054             meanalbedocld(j) = output_missing_value
1055             meanptop(j) = output_missing_value
1056             meantaucld(j) = output_missing_value
1057        end if
1058      enddo ! j
1059
1060      boxarea = 1./real(ncol)
1061     
1062      do 39 ibox=1,ncol
1063        do j=1,npoints
1064
1065          if (tau(j,ibox) .gt. (tauchk            )
1066     &      .and. ptop(j,ibox) .gt. 0.) then
1067              box_cloudy(j,ibox)=.true.
1068          endif
1069
1070          if (box_cloudy(j,ibox)) then
1071
1072              if (sunlit(j).eq.1 .or. top_height .eq. 3) then
1073
1074                boxtau(j,ibox) = tau(j,ibox)
1075
1076                if (tau(j,ibox) .ge. isccp_taumin) then
1077                   totalcldarea(j) = totalcldarea(j) + boxarea
1078               
1079                   !convert optical thickness to albedo
1080                   albedocld(j,ibox)
1081     &             = (tau(j,ibox)**0.895)/((tau(j,ibox)**0.895)+6.82)
1082         
1083                   !contribute to averaging
1084                   meanalbedocld(j) = meanalbedocld(j)
1085     &                                +albedocld(j,ibox)*boxarea
1086
1087                end if
1088
1089            endif
1090
1091          endif
1092
1093          if (sunlit(j).eq.1 .or. top_height .eq. 3) then
1094
1095           if (box_cloudy(j,ibox)) then
1096         
1097              !convert ptop to millibars
1098              ptop(j,ibox)=ptop(j,ibox) / 100.
1099           
1100              !save for output cloud top pressure and optical thickness
1101              boxptop(j,ibox) = ptop(j,ibox)
1102   
1103              if (tau(j,ibox) .ge. isccp_taumin) then
1104                meanptop(j) = meanptop(j) + ptop(j,ibox)*boxarea
1105              end if           
1106
1107              !reset itau(j), ipres(j)
1108              itau(j) = 0
1109              ipres(j) = 0
1110
1111              !determine optical depth category
1112              if (tau(j,ibox) .lt. isccp_taumin) then
1113                  itau(j)=1
1114              else if (tau(j,ibox) .ge. isccp_taumin
1115     &                                   
1116     &          .and. tau(j,ibox) .lt. 1.3) then
1117                itau(j)=2
1118              else if (tau(j,ibox) .ge. 1.3
1119     &          .and. tau(j,ibox) .lt. 3.6) then
1120                itau(j)=3
1121              else if (tau(j,ibox) .ge. 3.6
1122     &          .and. tau(j,ibox) .lt. 9.4) then
1123                  itau(j)=4
1124              else if (tau(j,ibox) .ge. 9.4
1125     &          .and. tau(j,ibox) .lt. 23.) then
1126                  itau(j)=5
1127              else if (tau(j,ibox) .ge. 23.
1128     &          .and. tau(j,ibox) .lt. 60.) then
1129                  itau(j)=6
1130              else if (tau(j,ibox) .ge. 60.) then
1131                  itau(j)=7
1132              end if
1133
1134              !determine cloud top pressure category
1135              if (    ptop(j,ibox) .gt. 0. 
1136     &          .and.ptop(j,ibox) .lt. 180.) then
1137                  ipres(j)=1
1138              else if(ptop(j,ibox) .ge. 180.
1139     &          .and.ptop(j,ibox) .lt. 310.) then
1140                  ipres(j)=2
1141              else if(ptop(j,ibox) .ge. 310.
1142     &          .and.ptop(j,ibox) .lt. 440.) then
1143                  ipres(j)=3
1144              else if(ptop(j,ibox) .ge. 440.
1145     &          .and.ptop(j,ibox) .lt. 560.) then
1146                  ipres(j)=4
1147              else if(ptop(j,ibox) .ge. 560.
1148     &          .and.ptop(j,ibox) .lt. 680.) then
1149                  ipres(j)=5
1150              else if(ptop(j,ibox) .ge. 680.
1151     &          .and.ptop(j,ibox) .lt. 800.) then
1152                  ipres(j)=6
1153              else if(ptop(j,ibox) .ge. 800.) then
1154                  ipres(j)=7
1155              end if
1156
1157              !update frequencies
1158              if(ipres(j) .gt. 0.and.itau(j) .gt. 0) then
1159              fq_isccp(j,itau(j),ipres(j))=
1160     &          fq_isccp(j,itau(j),ipres(j))+ boxarea
1161              end if
1162
1163            end if
1164
1165          end if
1166                       
1167        enddo ! j
116839    continue
1169     
1170      !compute mean cloud properties
1171      do j=1,npoints
1172        if (totalcldarea(j) .gt. 0.) then
1173          ! code above guarantees that totalcldarea > 0
1174          ! only if sunlit .eq. 1 .or. top_height = 3
1175          ! and applies only to clouds with tau > isccp_taumin
1176          meanptop(j) = meanptop(j) / totalcldarea(j)
1177          meanalbedocld(j) = meanalbedocld(j) / totalcldarea(j)
1178          meantaucld(j) = (6.82/((1./meanalbedocld(j))-1.))**(1./0.895)
1179        else
1180          ! this code is necessary so that in the case that totalcldarea = 0.,
1181          ! that these variables, which are in-cloud averages, are set to missing
1182          ! note that totalcldarea will be 0. if all the clouds in the grid box have
1183          ! tau < isccp_taumin
1184          meanptop(j) = output_missing_value
1185          meanalbedocld(j) = output_missing_value
1186          meantaucld(j) = output_missing_value
1187        end if
1188      enddo ! j
1189!
1190!     ---------------------------------------------------!
1191
1192!     ---------------------------------------------------!
1193!     OPTIONAL PRINTOUT OF DATA TO CHECK PROGRAM
1194!
1195      if (debugcol.ne.0) then
1196!     
1197         do j=1,npoints,debugcol
1198
1199            !produce character output
1200            do ilev=1,nlev
1201              do ibox=1,ncol
1202                   acc(ilev,ibox)=0
1203              enddo
1204            enddo
1205
1206            do ilev=1,nlev
1207              do ibox=1,ncol
1208                   acc(ilev,ibox)=frac_out(j,ibox,ilev)*2
1209                   if (levmatch(j,ibox) .eq. ilev)
1210     &                 acc(ilev,ibox)=acc(ilev,ibox)+1
1211              enddo
1212            enddo
1213
1214             !print test
1215
1216          write(ftn09,11) j
121711        format('ftn09.',i4.4)
1218          open(9, FILE=ftn09, FORM='FORMATTED')
1219
1220             write(9,'(a1)') ' '
1221             write(9,'(10i5)')
1222     &                  (ilev,ilev=5,nlev,5)
1223             write(9,'(a1)') ' '
1224             
1225             do ibox=1,ncol
1226               write(9,'(40(a1),1x,40(a1))')
1227     &           (cchar_realtops(acc(ilev,ibox)+1),ilev=1,nlev)
1228     &           ,(cchar(acc(ilev,ibox)+1),ilev=1,nlev)
1229             end do
1230             close(9)
1231
1232             if (ncolprint.ne.0) then
1233               write(6,'(a1)') ' '
1234                    write(6,'(a2,1X,5(a7,1X),a50)')
1235     &                  'ilev',
1236     &                  'pfull','at',
1237     &                  'cc*100','dem_s','dtau_s',
1238     &                  'cchar'
1239
1240!               do 4012 ilev=1,nlev
1241!                    write(6,'(60i2)') (box(i,ilev),i=1,ncolprint)
1242!                   write(6,'(i2,1X,5(f7.2,1X),50(a1))')
1243!     &                  ilev,
1244!     &                  pfull(j,ilev)/100.,at(j,ilev),
1245!     &                  cc(j,ilev)*100.0,dem_s(j,ilev),dtau_s(j,ilev)
1246!     &                  ,(cchar(acc(ilev,ibox)+1),ibox=1,ncolprint)
1247!4012           continue
1248               write (6,'(a)') 'skt(j):'
1249               write (6,'(8f7.2)') skt(j)
1250                                     
1251               write (6,'(8I7)') (ibox,ibox=1,ncolprint)
1252           
1253               write (6,'(a)') 'tau:'
1254               write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint)
1255   
1256               write (6,'(a)') 'tb:'
1257               write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint)
1258   
1259               write (6,'(a)') 'ptop:'
1260               write (6,'(8f7.2)') (ptop(j,ibox),ibox=1,ncolprint)
1261             endif
1262   
1263        enddo
1264       
1265      end if
1266
1267      return
1268      end
1269
1270
Note: See TracBrowser for help on using the repository browser.