source: LMDZ4/trunk/libf/cosp/icarus.F @ 2569

Last change on this file since 2569 was 1414, checked in by idelkadi, 14 years ago

Passage a la version cosp.v1.3 pour le Lidar et ISCCP
Corrections de bugs pour ISCCP et optimisation pour le Lidar

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