source: LMDZ4/branches/LMDZ4-dev-20091210/libf/cosp/icarus.F @ 1280

Last change on this file since 1280 was 1268, checked in by yann meurdesoif, 15 years ago

correction ligne commençant par !$ -> refuse de compiler en openMP

YM

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