Ignore:
Timestamp:
Mar 5, 2014, 2:19:12 PM (10 years ago)
Author:
lguez
Message:

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

File:
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/isccp_cloud_types.F90

    r1988 r1992  
    1 !
     1
    22! $Id $
    3 !
    4       SUBROUTINE ISCCP_CLOUD_TYPES(
    5      &     debug,
    6      &     debugcol,
    7      &     npoints,
    8      &     sunlit,
    9      &     nlev,
    10      &     ncol,
    11      &     seed,
    12      &     pfull,
    13      &     phalf,
    14      &     qv,
    15      &     cc,
    16      &     conv,
    17      &     dtau_s,
    18      &     dtau_c,
    19      &     top_height,
    20      &     overlap,
    21      &     tautab,
    22      &     invtau,
    23      &     skt,
    24      &     emsfc_lw,
    25      &     at,dem_s,dem_c,
    26      &     fq_isccp,
    27      &     totalcldarea,
    28      &     meanptop,
    29      &     meantaucld,
    30      &     boxtau,
    31      &     boxptop
    32      &)
    33        
    34 
    35 ! Copyright Steve Klein and Mark Webb 2002 - all rights reserved.
    36 !
    37 ! This code is available without charge with the following conditions:
    38 !
    39 !  1. The code is available for scientific purposes and is not for
    40 !     commercial use.
    41 !  2. Any improvements you make to the code should be made available
    42 !     to the to the authors for incorporation into a future release.
    43 !  3. The code should not be used in any way that brings the authors
    44 !     or their employers into disrepute.
    45 
    46       implicit none
    47 
    48 !     NOTE:   the maximum number of levels and columns is set by
    49 !             the following parameter statement
    50 
    51       INTEGER ncolprint
    52      
    53 !     -----
    54 !     Input
    55 !     -----
    56 
    57       INTEGER npoints                   !  number of model points in the horizontal
    58 c      PARAMETER(npoints=6722)
    59       INTEGER nlev                      !  number of model levels in column
    60       INTEGER ncol                      !  number of subcolumns
    61 
    62       INTEGER sunlit(npoints)           !  1 for day points, 0 for night time
    63 
    64       INTEGER seed(npoints)             !  seed value for random number generator
    65 c                                       !  ( see Numerical Recipes Chapter 7)
    66 c                                       !  It is recommended that the seed is set
    67 c                                       !  to a different value for each model
    68 c                                       !  gridbox it is called on, as it is
    69 c                                          !  possible that the choice of the samec
    70 c                                        !  seed value every time may introduce some
    71 c                                        !  statistical bias in the results, particularly
    72 c                                        !  for low values of NCOL.
    73 c
    74       REAL pfull(npoints,nlev)                      !  pressure of full model levels (Pascals)
    75 c                                        !  pfull(npoints,1)    is    top level of model
    76 c                                        !  pfull(npoints,nlev) is bottom level of model
    77 
    78       REAL phalf(npoints,nlev+1)        !  pressure of half model levels (Pascals)
    79 c                                        !  phalf(npoints,1)    is    top       of model
    80 c                                        !  phalf(npoints,nlev+1) is the surface pressure
    81 
    82       REAL qv(npoints,nlev)             !  water vapor specific humidity (kg vapor/ kg air)
    83 c                                        !         on full model levels
    84 
    85       REAL cc(npoints,nlev)             !  input cloud cover in each model level (fraction)
    86 c                                        !  NOTE:  This is the HORIZONTAL area of each
    87 c                                        !         grid box covered by clouds
    88 
    89       REAL conv(npoints,nlev)           !  input convective cloud cover in each model level (fraction)
    90 c                                        !  NOTE:  This is the HORIZONTAL area of each
    91 c                                        !         grid box covered by convective clouds
    92 
    93       REAL dtau_s(npoints,nlev)         !  mean 0.67 micron optical depth of stratiform
    94 c                                        !  clouds in each model level
    95 c                                        !  NOTE:  this the cloud optical depth of only the
    96 c                                        !         cloudy part of the grid box, it is not weighted
    97 c                                        !         with the 0 cloud optical depth of the clear
    98 c                                        !         part of the grid box
    99 
    100       REAL dtau_c(npoints,nlev)         !  mean 0.67 micron optical depth of convective
    101 c                                        !  clouds in each
    102 c                                        !  model level.  Same note applies as in dtau_s.
    103 
    104       INTEGER overlap                   !  overlap type
    105                                        
    106 !  1=max
    107                                        
    108 !  2=rand
    109 !  3=max/rand
    110 
    111       INTEGER top_height                !  1 = adjust top height using both a computed
    112 c                                        !  infrared brightness temperature and the visible
    113 c                                        !  optical depth to adjust cloud top pressure. Note
    114 c                                        !  that this calculation is most appropriate to compare
    115 c                                        !  to ISCCP data during sunlit hours.
    116 c                                        !  2 = do not adjust top height, that is cloud top
    117 c                                        !  pressure is the actual cloud top pressure
    118 c                                        !  in the model
    119 c                                        !  3 = adjust top height using only the computed
    120 c                                        !  infrared brightness temperature. Note that this
    121 c                                        !  calculation is most appropriate to compare to ISCCP
    122 c                                        !  IR only algortihm (i.e. you can compare to nighttime
    123 c                                        !  ISCCP data with this option)
    124 
    125       REAL tautab(0:255)                !  ISCCP table for converting count value to
    126 c                                        !  optical thickness
    127 
    128       INTEGER invtau(-20:45000)         !  ISCCP table for converting optical thickness
    129 c                                        !  to count value
    130 !
    131 !     The following input variables are used only if top_height = 1 or top_height = 3
    132 !
    133       REAL skt(npoints)                 !  skin Temperature (K)
    134       REAL emsfc_lw                     !  10.5 micron emissivity of surface (fraction)                                           
    135       REAL at(npoints,nlev)                   !  temperature in each model level (K)
    136       REAL dem_s(npoints,nlev)                !  10.5 micron longwave emissivity of stratiform
    137 c                                        !  clouds in each
    138 c                                        !  model level.  Same note applies as in dtau_s.
    139       REAL dem_c(npoints,nlev)                  !  10.5 micron longwave emissivity of convective
    140 c                                        !  clouds in each
    141 c                                        !  model level.  Same note applies as in dtau_s.
    142 cIM reg.dyn BEG
    143       REAL t1, t2
    144 !     REAL w(npoints)                   !vertical wind at 500 hPa
    145 !     LOGICAL pct_ocean(npoints)        !TRUE if oceanic point, FALSE otherway
    146 !     INTEGER iw(npoints) , nw
    147 !     REAL wmin, pas_w
    148 !     INTEGER k, l, iwmx
    149 !     PARAMETER(wmin=-100.,pas_w=10.,iwmx=30)
    150 !     REAL fq_dynreg(7,7,iwmx)
    151 !     REAL nfq_dynreg(7,7,iwmx)
    152 !     LOGICAL pctj(7,7,iwmx)
    153 cIM reg.dyn END
    154 !     ------
    155 !     Output
    156 !     ------
    157 
    158       REAL fq_isccp(npoints,7,7)        !  the fraction of the model grid box covered by
    159 c                                        !  each of the 49 ISCCP D level cloud types
    160 
    161       REAL totalcldarea(npoints)        !  the fraction of model grid box columns
    162 c                                        !  with cloud somewhere in them.  This should
    163 c                                        !  equal the sum over all entries of fq_isccp
    164        
    165        
    166 c      ! The following three means are averages over the cloudy areas only.  If no
    167 c      ! clouds are in grid box all three quantities should equal zero.       
    168                                        
    169       REAL meanptop(npoints)            !  mean cloud top pressure (mb) - linear averaging
    170 c                                        !  in cloud top pressure.
    171                                        
    172       REAL meantaucld(npoints)          !  mean optical thickness
    173 c                                        !  linear averaging in albedo performed.
    174      
    175       REAL boxtau(npoints,ncol)         !  optical thickness in each column
    176      
    177       REAL boxptop(npoints,ncol)        !  cloud top pressure (mb) in each column
    178                                        
    179                                                                                                                        
    180 !
    181 !     ------
    182 !     Working variables added when program updated to mimic Mark Webb's PV-Wave code
    183 !     ------
    184 
    185       REAL frac_out(npoints,ncol,nlev) ! boxes gridbox divided up into
    186 c                                        ! Equivalent of BOX in original version, but
    187 c                                        ! indexed by column then row, rather than
    188 c                                        ! by row then column
    189 
    190       REAL tca(npoints,0:nlev) ! total cloud cover in each model level (fraction)
    191 c                                        ! with extra layer of zeroes on top
    192 c                                        ! in this version this just contains the values input
    193 c                                        ! from cc but with an extra level
    194       REAL cca(npoints,nlev)         ! convective cloud cover in each model level (fraction)
    195 c                                        ! from conv
    196 
    197       REAL threshold(npoints,ncol)   ! pointer to position in gridbox
    198       REAL maxocc(npoints,ncol)      ! Flag for max overlapped conv cld
    199       REAL maxosc(npoints,ncol)      ! Flag for max overlapped strat cld
    200 
    201       REAL boxpos(npoints,ncol)      ! ordered pointer to position in gridbox
    202 
    203       REAL threshold_min(npoints,ncol) ! minimum value to define range in with new threshold
    204 c                                        ! is chosen
    205 
    206       REAL dem(npoints,ncol),bb(npoints)     !  working variables for 10.5 micron longwave
    207 c                                        !  emissivity in part of
    208 c                                        !  gridbox under consideration
    209 
    210       REAL ran(npoints)                   ! vector of random numbers
    211       REAL ptrop(npoints)
    212       REAL attrop(npoints)
    213       REAL attropmin (npoints)
    214       REAL atmax(npoints)
    215       REAL atmin(npoints)
    216       REAL btcmin(npoints)
    217       REAL transmax(npoints)
    218 
    219       INTEGER i,j,ilev,ibox,itrop(npoints)
    220       INTEGER ipres(npoints)
    221       INTEGER itau(npoints),ilev2
    222       INTEGER acc(nlev,ncol)
    223       INTEGER match(npoints,nlev-1)
    224       INTEGER nmatch(npoints)
    225       INTEGER levmatch(npoints,ncol)
    226      
    227 c      !variables needed for water vapor continuum absorption
    228       real fluxtop_clrsky(npoints),trans_layers_above_clrsky(npoints)
    229       real taumin(npoints)
    230       real dem_wv(npoints,nlev), wtmair, wtmh20, Navo, grav, pstd, t0
    231       real press(npoints), dpress(npoints), atmden(npoints)
    232       real rvh20(npoints), wk(npoints), rhoave(npoints)
    233       real rh20s(npoints), rfrgn(npoints)
    234       real tmpexp(npoints),tauwv(npoints)
    235      
    236       character*1 cchar(6),cchar_realtops(6)
    237       integer icycle
    238       REAL tau(npoints,ncol)
    239       LOGICAL box_cloudy(npoints,ncol)
    240       REAL tb(npoints,ncol)
    241       REAL ptop(npoints,ncol)
    242       REAL emcld(npoints,ncol)
    243       REAL fluxtop(npoints,ncol)
    244       REAL trans_layers_above(npoints,ncol)
    245       real isccp_taumin,fluxtopinit(npoints),tauir(npoints)
    246       real meanalbedocld(npoints)
    247       REAL albedocld(npoints,ncol)
    248       real boxarea
    249       integer debug       ! set to non-zero value to print out inputs
    250 c                          ! with step debug
    251       integer debugcol    ! set to non-zero value to print out column
    252 c                          ! decomposition with step debugcol
    253 
    254       integer index1(npoints),num1,jj
    255       real rec2p13,tauchk
    256 
    257       character*10 ftn09
    258      
    259       DATA isccp_taumin / 0.3 /
    260       DATA cchar / ' ','-','1','+','I','+'/
    261       DATA cchar_realtops / ' ',' ','1','1','I','I'/
    262 
    263       tauchk = -1.*log(0.9999999)
    264       rec2p13=1./2.13
    265 
    266       ncolprint=0
    267 
    268 cIM
    269 c     PRINT*,' isccp_cloud_types npoints=',npoints
    270 c
    271 c      if ( debug.ne.0 ) then
    272 c          j=1
    273 c          write(6,'(a10)') 'j='
    274 c          write(6,'(8I10)') j
    275 c          write(6,'(a10)') 'debug='
    276 c          write(6,'(8I10)') debug
    277 c          write(6,'(a10)') 'debugcol='
    278 c          write(6,'(8I10)') debugcol
    279 c          write(6,'(a10)') 'npoints='
    280 c          write(6,'(8I10)') npoints
    281 c          write(6,'(a10)') 'nlev='
    282 c          write(6,'(8I10)') nlev
    283 c          write(6,'(a10)') 'ncol='
    284 c          write(6,'(8I10)') ncol
    285 c          write(6,'(a10)') 'top_height='
    286 c          write(6,'(8I10)') top_height
    287 c          write(6,'(a10)') 'overlap='
    288 c          write(6,'(8I10)') overlap
    289 c          write(6,'(a10)') 'emsfc_lw='
    290 c          write(6,'(8f10.2)') emsfc_lw
    291 c          write(6,'(a10)') 'tautab='
    292 c          write(6,'(8f10.2)') tautab
    293 c          write(6,'(a10)') 'invtau(1:100)='
    294 c          write(6,'(8i10)') (invtau(i),i=1,100)
    295 c          do j=1,npoints,debug
    296 c          write(6,'(a10)') 'j='
    297 c          write(6,'(8I10)') j
    298 c          write(6,'(a10)') 'sunlit='
    299 c          write(6,'(8I10)') sunlit(j)
    300 c          write(6,'(a10)') 'seed='
    301 c          write(6,'(8I10)') seed(j)
    302 c          write(6,'(a10)') 'pfull='
    303 c          write(6,'(8f10.2)') (pfull(j,i),i=1,nlev)
    304 c          write(6,'(a10)') 'phalf='
    305 c          write(6,'(8f10.2)') (phalf(j,i),i=1,nlev+1)
    306 c          write(6,'(a10)') 'qv='
    307 c          write(6,'(8f10.3)') (qv(j,i),i=1,nlev)
    308 c          write(6,'(a10)') 'cc='
    309 c          write(6,'(8f10.3)') (cc(j,i),i=1,nlev)
    310 c          write(6,'(a10)') 'conv='
    311 c          write(6,'(8f10.2)') (conv(j,i),i=1,nlev)
    312 c          write(6,'(a10)') 'dtau_s='
    313 c          write(6,'(8g12.5)') (dtau_s(j,i),i=1,nlev)
    314 c          write(6,'(a10)') 'dtau_c='
    315 c          write(6,'(8f10.2)') (dtau_c(j,i),i=1,nlev)
    316 c          write(6,'(a10)') 'skt='
    317 c          write(6,'(8f10.2)') skt(j)
    318 c          write(6,'(a10)') 'at='
    319 c          write(6,'(8f10.2)') (at(j,i),i=1,nlev)
    320 c          write(6,'(a10)') 'dem_s='
    321 c          write(6,'(8f10.3)') (dem_s(j,i),i=1,nlev)
    322 c          write(6,'(a10)') 'dem_c='
    323 c          write(6,'(8f10.2)') (dem_c(j,i),i=1,nlev)
    324 c          enddo
    325 c      endif
    326 
    327 !     ---------------------------------------------------!
    328 
    329 !     assign 2d tca array using 1d input array cc
    330 
    331       do j=1,npoints
    332         tca(j,0)=0
    333       enddo
    334  
    335       do ilev=1,nlev
    336         do j=1,npoints
    337           tca(j,ilev)=cc(j,ilev)
    338         enddo
    339       enddo
    340 
    341 !     assign 2d cca array using 1d input array conv
    342 
    343       do ilev=1,nlev
    344 cIM pas besoin        do ibox=1,ncol
    345           do j=1,npoints
    346             cca(j,ilev)=conv(j,ilev)
    347           enddo
    348 cIM        enddo
    349       enddo
    350 
    351 cIM
    352 !     do j=1, iwmx
    353 !     do l=1, 7
    354 !     do k=1, 7
    355 !       fq_dynreg(k,l,j) =0.
    356 !       nfq_dynreg(k,l,j) =0.
    357 !      enddo !k
    358 !     enddo !l
    359 !     enddo !j
    360 cIM
    361 cIM
    362 c      if (ncolprint.ne.0) then
    363 c        do j=1,npoints,1000
    364 c        write(6,'(a10)') 'j='
    365 c        write(6,'(8I10)') j
    366 c        write (6,'(a)') 'seed:'
    367 c        write (6,'(I3.2)') seed(j)
    368 
    369 c        write (6,'(a)') 'tca_pp_rev:'
    370 c        write (6,'(8f5.2)')
    371 c     &   ((tca(j,ilev)),
    372 c     &      ilev=1,nlev)
    373 
    374 c        write (6,'(a)') 'cca_pp_rev:'
    375 c        write (6,'(8f5.2)')
    376 c     &   ((cca(j,ilev),ibox=1,ncolprint),ilev=1,nlev)
    377 c        enddo
    378 c      endif
    379 
    380       if (top_height .eq. 1 .or. top_height .eq. 3) then
    381 
    382       do j=1,npoints
    383           ptrop(j)=5000.
    384           atmin(j) = 400.
    385           attropmin(j) = 400.
    386           atmax(j) = 0.
    387           attrop(j) = 120.
    388           itrop(j) = 1
    389       enddo
    390 
    391       do 12 ilev=1,nlev
    392         do j=1,npoints
    393          if (pfull(j,ilev) .lt. 40000. .and.
    394      &          pfull(j,ilev) .gt.  5000. .and.
    395      &          at(j,ilev) .lt. attropmin(j)) then
    396                 ptrop(j) = pfull(j,ilev)
    397                 attropmin(j) = at(j,ilev)
    398                 attrop(j) = attropmin(j)
    399                 itrop(j)=ilev
    400            end if
    401            if (at(j,ilev) .gt. atmax(j)) atmax(j)=at(j,ilev)
    402            if (at(j,ilev) .lt. atmin(j)) atmin(j)=at(j,ilev)
    403         enddo
    404 12    continue
    405 
    406       end if
    407 
    408 !     -----------------------------------------------------!
    409 
    410 !     ---------------------------------------------------!
    411 
    412 cIM
    413 c     do 13 ilev=1,nlev
    414 cnum1=0
    415 c       do j=1,npoints
    416 c           if (cc(j,ilev) .lt. 0. .or. cc(j,ilev) .gt. 1.) then
    417 c        num1=num1+1
    418 c        index1(num1)=j
    419 c           end if
    420 c       enddo
    421 c       do jj=1,num1   
    422 c        j=index1(jj)
    423 c               write(6,*)  ' error = cloud fraction less than zero'
    424 c        write(6,*) ' or '
    425 c               write(6,*)  ' error = cloud fraction greater than 1'
    426 c        write(6,*) 'value at point ',j,' is ',cc(j,ilev)
    427 c        write(6,*) 'level ',ilev
    428 c               STOP
    429 c       enddo
    430 cnum1=0
    431 c       do j=1,npoints
    432 c           if (conv(j,ilev) .lt. 0. .or. conv(j,ilev) .gt. 1.) then
    433 c        num1=num1+1
    434 c        index1(num1)=j
    435 c           end if
    436 c       enddo
    437 c       do jj=1,num1   
    438 c        j=index1(jj)
    439 c               write(6,*) 
    440 c    &           ' error = convective cloud fraction less than zero'
    441 c        write(6,*) ' or '
    442 c               write(6,*) 
    443 c    &           ' error = convective cloud fraction greater than 1'
    444 c        write(6,*) 'value at point ',j,' is ',conv(j,ilev)
    445 c        write(6,*) 'level ',ilev
    446 c               STOP
    447 c       enddo
    448 
    449 cnum1=0
    450 c       do j=1,npoints
    451 c           if (dtau_s(j,ilev) .lt. 0.) then
    452 c        num1=num1+1
    453 c        index1(num1)=j
    454 c           end if
    455 c       enddo
    456 c       do jj=1,num1   
    457 c        j=index1(jj)
    458 c               write(6,*) 
    459 c    &           ' error = stratiform cloud opt. depth less than zero'
    460 c        write(6,*) 'value at point ',j,' is ',dtau_s(j,ilev)
    461 c        write(6,*) 'level ',ilev
    462 c               STOP
    463 c       enddo
    464 cnum1=0
    465 c       do j=1,npoints
    466 c           if (dtau_c(j,ilev) .lt. 0.) then
    467 c        num1=num1+1
    468 c        index1(num1)=j
    469 c           end if
    470 c       enddo
    471 c       do jj=1,num1   
    472 c        j=index1(jj)
    473 c               write(6,*) 
    474 c    &           ' error = convective cloud opt. depth less than zero'
    475 c        write(6,*) 'value at point ',j,' is ',dtau_c(j,ilev)
    476 c        write(6,*) 'level ',ilev
    477 c               STOP
    478 c       enddo
    479 
    480 cnum1=0
    481 c       do j=1,npoints
    482 c           if (dem_s(j,ilev) .lt. 0. .or. dem_s(j,ilev) .gt. 1.) then
    483 c        num1=num1+1
    484 c        index1(num1)=j
    485 c           end if
    486 c       enddo
    487 c       do jj=1,num1   
    488 c        j=index1(jj)
    489 c               write(6,*) 
    490 c    &           ' error = stratiform cloud emissivity less than zero'
    491 c        write(6,*)'or'
    492 c               write(6,*) 
    493 c    &           ' error = stratiform cloud emissivity greater than 1'
    494 c        write(6,*) 'value at point ',j,' is ',dem_s(j,ilev)
    495 c        write(6,*) 'level ',ilev
    496 c               STOP
    497 c       enddo
    498 
    499 cnum1=0
    500 c       do j=1,npoints
    501 c           if (dem_c(j,ilev) .lt. 0. .or. dem_c(j,ilev) .gt. 1.) then
    502 c        num1=num1+1
    503 c        index1(num1)=j
    504 c           end if
    505 c       enddo
    506 c       do jj=1,num1   
    507 c        j=index1(jj)
    508 c               write(6,*) 
    509 c    &           ' error = convective cloud emissivity less than zero'
    510 c        write(6,*)'or'
    511 c               write(6,*) 
    512 c    &           ' error = convective cloud emissivity greater than 1'
    513 c               write (6,*)
    514 c    &          'j=',j,'ilev=',ilev,'dem_c(j,ilev) =',dem_c(j,ilev)
    515 c               STOP
    516 c       enddo
    517 c13    continue
    518 
    519 
    520       do ibox=1,ncol
    521         do j=1,npoints
    522           boxpos(j,ibox)=(ibox-.5)/ncol
    523         enddo
    524       enddo
    525 
    526 !     ---------------------------------------------------!
    527 !     Initialise working variables
    528 !     ---------------------------------------------------!
    529 
    530 !     Initialised frac_out to zero
    531 
    532       do ilev=1,nlev
    533         do ibox=1,ncol
    534           do j=1,npoints
    535             frac_out(j,ibox,ilev)=0.0
    536           enddo
    537         enddo
    538       enddo
    539 
    540 cIM
    541 c      if (ncolprint.ne.0) then
    542 c        write (6,'(a)') 'frac_out_pp_rev:'
    543 c          do j=1,npoints,1000
    544 c          write(6,'(a10)') 'j='
    545 c          write(6,'(8I10)') j
    546 c          write (6,'(8f5.2)')
    547 c     &     ((frac_out(j,ibox,ilev),ibox=1,ncolprint),ilev=1,nlev)
    548 
    549 c          enddo
    550 c        write (6,'(a)') 'ncol:'
    551 c        write (6,'(I3)') ncol
    552 c      endif
    553 c      if (ncolprint.ne.0) then
    554 c        write (6,'(a)') 'last_frac_pp:'
    555 c          do j=1,npoints,1000
    556 c          write(6,'(a10)') 'j='
    557 c          write(6,'(8I10)') j
    558 c          write (6,'(8f5.2)') (tca(j,0))
    559 c          enddo
    560 c      endif
    561 
    562 !     ---------------------------------------------------!
    563 !     ALLOCATE CLOUD INTO BOXES, FOR NCOLUMNS, NLEVELS
    564 !     frac_out is the array that contains the information
    565 !     where 0 is no cloud, 1 is a stratiform cloud and 2 is a
    566 !     convective cloud
    567      
    568       !loop over vertical levels
    569       DO 200 ilev = 1,nlev
    570                                  
    571 !     Initialise threshold
    572 
    573         IF (ilev.eq.1) then
    574             ! If max overlap
    575             IF (overlap.eq.1) then
    576               ! select pixels spread evenly
    577               ! across the gridbox
    578               DO ibox=1,ncol
    579                 do j=1,npoints
    580                   threshold(j,ibox)=boxpos(j,ibox)
    581                 enddo
    582               enddo
     3
     4SUBROUTINE isccp_cloud_types(debug, debugcol, npoints, sunlit, nlev, ncol, &
     5    seed, pfull, phalf, qv, cc, conv, dtau_s, dtau_c, top_height, overlap, &
     6    tautab, invtau, skt, emsfc_lw, at, dem_s, dem_c, fq_isccp, totalcldarea, &
     7    meanptop, meantaucld, boxtau, boxptop)
     8
     9
     10  ! Copyright Steve Klein and Mark Webb 2002 - all rights reserved.
     11
     12  ! This code is available without charge with the following conditions:
     13
     14  ! 1. The code is available for scientific purposes and is not for
     15  ! commercial use.
     16  ! 2. Any improvements you make to the code should be made available
     17  ! to the to the authors for incorporation into a future release.
     18  ! 3. The code should not be used in any way that brings the authors
     19  ! or their employers into disrepute.
     20
     21  IMPLICIT NONE
     22
     23  ! NOTE:   the maximum number of levels and columns is set by
     24  ! the following parameter statement
     25
     26  INTEGER ncolprint
     27
     28  ! -----
     29  ! Input
     30  ! -----
     31
     32  INTEGER npoints !  number of model points in the horizontal
     33  ! PARAMETER(npoints=6722)
     34  INTEGER nlev !  number of model levels in column
     35  INTEGER ncol !  number of subcolumns
     36
     37  INTEGER sunlit(npoints) !  1 for day points, 0 for night time
     38
     39  INTEGER seed(npoints) !  seed value for random number generator
     40  ! !  ( see Numerical Recipes Chapter 7)
     41  ! !  It is recommended that the seed is set
     42  ! !  to a different value for each model
     43  ! !  gridbox it is called on, as it is
     44  ! !  possible that the choice of the samec
     45  ! !  seed value every time may introduce some
     46  ! !  statistical bias in the results, particularly
     47  ! !  for low values of NCOL.
     48
     49  REAL pfull(npoints, nlev) !  pressure of full model levels (Pascals)
     50  ! !  pfull(npoints,1)    is    top level of model
     51  ! !  pfull(npoints,nlev) is bottom level of model
     52
     53  REAL phalf(npoints, nlev+1) !  pressure of half model levels (Pascals)
     54  ! !  phalf(npoints,1)    is    top       of model
     55  ! !  phalf(npoints,nlev+1) is the surface pressure
     56
     57  REAL qv(npoints, nlev) !  water vapor specific humidity (kg vapor/ kg air)
     58  ! !         on full model levels
     59
     60  REAL cc(npoints, nlev) !  input cloud cover in each model level (fraction)
     61  ! !  NOTE:  This is the HORIZONTAL area of each
     62  ! !         grid box covered by clouds
     63
     64  REAL conv(npoints, nlev) !  input convective cloud cover in each model level (fraction)
     65  ! !  NOTE:  This is the HORIZONTAL area of each
     66  ! !         grid box covered by convective clouds
     67
     68  REAL dtau_s(npoints, nlev) !  mean 0.67 micron optical depth of stratiform
     69  ! !  clouds in each model level
     70  ! !  NOTE:  this the cloud optical depth of only the
     71  ! !         cloudy part of the grid box, it is not weighted
     72  ! !         with the 0 cloud optical depth of the clear
     73  ! !         part of the grid box
     74
     75  REAL dtau_c(npoints, nlev) !  mean 0.67 micron optical depth of convective
     76  ! !  clouds in each
     77  ! !  model level.  Same note applies as in dtau_s.
     78
     79  INTEGER overlap !  overlap type
     80
     81  ! 1=max
     82
     83  ! 2=rand
     84  ! 3=max/rand
     85
     86  INTEGER top_height !  1 = adjust top height using both a computed
     87  ! !  infrared brightness temperature and the visible
     88  ! !  optical depth to adjust cloud top pressure. Note
     89  ! !  that this calculation is most appropriate to compare
     90  ! !  to ISCCP data during sunlit hours.
     91  ! !  2 = do not adjust top height, that is cloud top
     92  ! !  pressure is the actual cloud top pressure
     93  ! !  in the model
     94  ! !  3 = adjust top height using only the computed
     95  ! !  infrared brightness temperature. Note that this
     96  ! !  calculation is most appropriate to compare to ISCCP
     97  ! !  IR only algortihm (i.e. you can compare to nighttime
     98  ! !  ISCCP data with this option)
     99
     100  REAL tautab(0:255) !  ISCCP table for converting count value to
     101  ! !  optical thickness
     102
     103  INTEGER invtau(-20:45000) !  ISCCP table for converting optical thickness
     104  ! !  to count value
     105
     106  ! The following input variables are used only if top_height = 1 or
     107  ! top_height = 3
     108
     109  REAL skt(npoints) !  skin Temperature (K)
     110  REAL emsfc_lw !  10.5 micron emissivity of surface (fraction)
     111  REAL at(npoints, nlev) !  temperature in each model level (K)
     112  REAL dem_s(npoints, nlev) !  10.5 micron longwave emissivity of stratiform
     113  ! !  clouds in each
     114  ! !  model level.  Same note applies as in dtau_s.
     115  REAL dem_c(npoints, nlev) !  10.5 micron longwave emissivity of convective
     116  ! !  clouds in each
     117  ! !  model level.  Same note applies as in dtau_s.
     118  ! IM reg.dyn BEG
     119  REAL t1, t2
     120  ! REAL w(npoints)                   !vertical wind at 500 hPa
     121  ! LOGICAL pct_ocean(npoints)        !TRUE if oceanic point, FALSE otherway
     122  ! INTEGER iw(npoints) , nw
     123  ! REAL wmin, pas_w
     124  ! INTEGER k, l, iwmx
     125  ! PARAMETER(wmin=-100.,pas_w=10.,iwmx=30)
     126  ! REAL fq_dynreg(7,7,iwmx)
     127  ! REAL nfq_dynreg(7,7,iwmx)
     128  ! LOGICAL pctj(7,7,iwmx)
     129  ! IM reg.dyn END
     130  ! ------
     131  ! Output
     132  ! ------
     133
     134  REAL fq_isccp(npoints, 7, 7) !  the fraction of the model grid box covered by
     135  ! !  each of the 49 ISCCP D level cloud types
     136
     137  REAL totalcldarea(npoints) !  the fraction of model grid box columns
     138  ! !  with cloud somewhere in them.  This should
     139  ! !  equal the sum over all entries of fq_isccp
     140
     141
     142  ! ! The following three means are averages over the cloudy areas only.  If
     143  ! no
     144  ! ! clouds are in grid box all three quantities should equal zero.
     145
     146  REAL meanptop(npoints) !  mean cloud top pressure (mb) - linear averaging
     147  ! !  in cloud top pressure.
     148
     149  REAL meantaucld(npoints) !  mean optical thickness
     150  ! !  linear averaging in albedo performed.
     151
     152  REAL boxtau(npoints, ncol) !  optical thickness in each column
     153
     154  REAL boxptop(npoints, ncol) !  cloud top pressure (mb) in each column
     155
     156
     157
     158  ! ------
     159  ! Working variables added when program updated to mimic Mark Webb's PV-Wave
     160  ! code
     161  ! ------
     162
     163  REAL frac_out(npoints, ncol, nlev) ! boxes gridbox divided up into
     164  ! ! Equivalent of BOX in original version, but
     165  ! ! indexed by column then row, rather than
     166  ! ! by row then column
     167
     168  REAL tca(npoints, 0:nlev) ! total cloud cover in each model level (fraction)
     169  ! ! with extra layer of zeroes on top
     170  ! ! in this version this just contains the values input
     171  ! ! from cc but with an extra level
     172  REAL cca(npoints, nlev) ! convective cloud cover in each model level (fraction)
     173  ! ! from conv
     174
     175  REAL threshold(npoints, ncol) ! pointer to position in gridbox
     176  REAL maxocc(npoints, ncol) ! Flag for max overlapped conv cld
     177  REAL maxosc(npoints, ncol) ! Flag for max overlapped strat cld
     178
     179  REAL boxpos(npoints, ncol) ! ordered pointer to position in gridbox
     180
     181  REAL threshold_min(npoints, ncol) ! minimum value to define range in with new threshold
     182  ! ! is chosen
     183
     184  REAL dem(npoints, ncol), bb(npoints) !  working variables for 10.5 micron longwave
     185  ! !  emissivity in part of
     186  ! !  gridbox under consideration
     187
     188  REAL ran(npoints) ! vector of random numbers
     189  REAL ptrop(npoints)
     190  REAL attrop(npoints)
     191  REAL attropmin(npoints)
     192  REAL atmax(npoints)
     193  REAL atmin(npoints)
     194  REAL btcmin(npoints)
     195  REAL transmax(npoints)
     196
     197  INTEGER i, j, ilev, ibox, itrop(npoints)
     198  INTEGER ipres(npoints)
     199  INTEGER itau(npoints), ilev2
     200  INTEGER acc(nlev, ncol)
     201  INTEGER match(npoints, nlev-1)
     202  INTEGER nmatch(npoints)
     203  INTEGER levmatch(npoints, ncol)
     204
     205  ! !variables needed for water vapor continuum absorption
     206  REAL fluxtop_clrsky(npoints), trans_layers_above_clrsky(npoints)
     207  REAL taumin(npoints)
     208  REAL dem_wv(npoints, nlev), wtmair, wtmh20, navo, grav, pstd, t0
     209  REAL press(npoints), dpress(npoints), atmden(npoints)
     210  REAL rvh20(npoints), wk(npoints), rhoave(npoints)
     211  REAL rh20s(npoints), rfrgn(npoints)
     212  REAL tmpexp(npoints), tauwv(npoints)
     213
     214  CHARACTER *1 cchar(6), cchar_realtops(6)
     215  INTEGER icycle
     216  REAL tau(npoints, ncol)
     217  LOGICAL box_cloudy(npoints, ncol)
     218  REAL tb(npoints, ncol)
     219  REAL ptop(npoints, ncol)
     220  REAL emcld(npoints, ncol)
     221  REAL fluxtop(npoints, ncol)
     222  REAL trans_layers_above(npoints, ncol)
     223  REAL isccp_taumin, fluxtopinit(npoints), tauir(npoints)
     224  REAL meanalbedocld(npoints)
     225  REAL albedocld(npoints, ncol)
     226  REAL boxarea
     227  INTEGER debug ! set to non-zero value to print out inputs
     228  ! ! with step debug
     229  INTEGER debugcol ! set to non-zero value to print out column
     230  ! ! decomposition with step debugcol
     231
     232  INTEGER index1(npoints), num1, jj
     233  REAL rec2p13, tauchk
     234
     235  CHARACTER *10 ftn09
     236
     237  DATA isccp_taumin/0.3/
     238  DATA cchar/' ', '-', '1', '+', 'I', '+'/
     239  DATA cchar_realtops/' ', ' ', '1', '1', 'I', 'I'/
     240
     241  tauchk = -1.*log(0.9999999)
     242  rec2p13 = 1./2.13
     243
     244  ncolprint = 0
     245
     246  ! IM
     247  ! PRINT*,' isccp_cloud_types npoints=',npoints
     248
     249  ! if ( debug.ne.0 ) then
     250  ! j=1
     251  ! write(6,'(a10)') 'j='
     252  ! write(6,'(8I10)') j
     253  ! write(6,'(a10)') 'debug='
     254  ! write(6,'(8I10)') debug
     255  ! write(6,'(a10)') 'debugcol='
     256  ! write(6,'(8I10)') debugcol
     257  ! write(6,'(a10)') 'npoints='
     258  ! write(6,'(8I10)') npoints
     259  ! write(6,'(a10)') 'nlev='
     260  ! write(6,'(8I10)') nlev
     261  ! write(6,'(a10)') 'ncol='
     262  ! write(6,'(8I10)') ncol
     263  ! write(6,'(a10)') 'top_height='
     264  ! write(6,'(8I10)') top_height
     265  ! write(6,'(a10)') 'overlap='
     266  ! write(6,'(8I10)') overlap
     267  ! write(6,'(a10)') 'emsfc_lw='
     268  ! write(6,'(8f10.2)') emsfc_lw
     269  ! write(6,'(a10)') 'tautab='
     270  ! write(6,'(8f10.2)') tautab
     271  ! write(6,'(a10)') 'invtau(1:100)='
     272  ! write(6,'(8i10)') (invtau(i),i=1,100)
     273  ! do j=1,npoints,debug
     274  ! write(6,'(a10)') 'j='
     275  ! write(6,'(8I10)') j
     276  ! write(6,'(a10)') 'sunlit='
     277  ! write(6,'(8I10)') sunlit(j)
     278  ! write(6,'(a10)') 'seed='
     279  ! write(6,'(8I10)') seed(j)
     280  ! write(6,'(a10)') 'pfull='
     281  ! write(6,'(8f10.2)') (pfull(j,i),i=1,nlev)
     282  ! write(6,'(a10)') 'phalf='
     283  ! write(6,'(8f10.2)') (phalf(j,i),i=1,nlev+1)
     284  ! write(6,'(a10)') 'qv='
     285  ! write(6,'(8f10.3)') (qv(j,i),i=1,nlev)
     286  ! write(6,'(a10)') 'cc='
     287  ! write(6,'(8f10.3)') (cc(j,i),i=1,nlev)
     288  ! write(6,'(a10)') 'conv='
     289  ! write(6,'(8f10.2)') (conv(j,i),i=1,nlev)
     290  ! write(6,'(a10)') 'dtau_s='
     291  ! write(6,'(8g12.5)') (dtau_s(j,i),i=1,nlev)
     292  ! write(6,'(a10)') 'dtau_c='
     293  ! write(6,'(8f10.2)') (dtau_c(j,i),i=1,nlev)
     294  ! write(6,'(a10)') 'skt='
     295  ! write(6,'(8f10.2)') skt(j)
     296  ! write(6,'(a10)') 'at='
     297  ! write(6,'(8f10.2)') (at(j,i),i=1,nlev)
     298  ! write(6,'(a10)') 'dem_s='
     299  ! write(6,'(8f10.3)') (dem_s(j,i),i=1,nlev)
     300  ! write(6,'(a10)') 'dem_c='
     301  ! write(6,'(8f10.2)') (dem_c(j,i),i=1,nlev)
     302  ! enddo
     303  ! endif
     304
     305  ! ---------------------------------------------------!
     306
     307  ! assign 2d tca array using 1d input array cc
     308
     309  DO j = 1, npoints
     310    tca(j, 0) = 0
     311  END DO
     312
     313  DO ilev = 1, nlev
     314    DO j = 1, npoints
     315      tca(j, ilev) = cc(j, ilev)
     316    END DO
     317  END DO
     318
     319  ! assign 2d cca array using 1d input array conv
     320
     321  DO ilev = 1, nlev
     322    ! IM pas besoin        do ibox=1,ncol
     323    DO j = 1, npoints
     324      cca(j, ilev) = conv(j, ilev)
     325    END DO
     326    ! IM        enddo
     327  END DO
     328
     329  ! IM
     330  ! do j=1, iwmx
     331  ! do l=1, 7
     332  ! do k=1, 7
     333  ! fq_dynreg(k,l,j) =0.
     334  ! nfq_dynreg(k,l,j) =0.
     335  ! enddo !k
     336  ! enddo !l
     337  ! enddo !j
     338  ! IM
     339  ! IM
     340  ! if (ncolprint.ne.0) then
     341  ! do j=1,npoints,1000
     342  ! write(6,'(a10)') 'j='
     343  ! write(6,'(8I10)') j
     344  ! write (6,'(a)') 'seed:'
     345  ! write (6,'(I3.2)') seed(j)
     346
     347  ! write (6,'(a)') 'tca_pp_rev:'
     348  ! write (6,'(8f5.2)')
     349  ! &   ((tca(j,ilev)),
     350  ! &      ilev=1,nlev)
     351
     352  ! write (6,'(a)') 'cca_pp_rev:'
     353  ! write (6,'(8f5.2)')
     354  ! &   ((cca(j,ilev),ibox=1,ncolprint),ilev=1,nlev)
     355  ! enddo
     356  ! endif
     357
     358  IF (top_height==1 .OR. top_height==3) THEN
     359
     360    DO j = 1, npoints
     361      ptrop(j) = 5000.
     362      atmin(j) = 400.
     363      attropmin(j) = 400.
     364      atmax(j) = 0.
     365      attrop(j) = 120.
     366      itrop(j) = 1
     367    END DO
     368
     369    DO ilev = 1, nlev
     370      DO j = 1, npoints
     371        IF (pfull(j,ilev)<40000. .AND. pfull(j,ilev)>5000. .AND. &
     372            at(j,ilev)<attropmin(j)) THEN
     373          ptrop(j) = pfull(j, ilev)
     374          attropmin(j) = at(j, ilev)
     375          attrop(j) = attropmin(j)
     376          itrop(j) = ilev
     377        END IF
     378        IF (at(j,ilev)>atmax(j)) atmax(j) = at(j, ilev)
     379        IF (at(j,ilev)<atmin(j)) atmin(j) = at(j, ilev)
     380      END DO
     381    END DO
     382
     383  END IF
     384
     385  ! -----------------------------------------------------!
     386
     387  ! ---------------------------------------------------!
     388
     389  ! IM
     390  ! do 13 ilev=1,nlev
     391  ! num1=0
     392  ! do j=1,npoints
     393  ! if (cc(j,ilev) .lt. 0. .or. cc(j,ilev) .gt. 1.) then
     394  ! num1=num1+1
     395  ! index1(num1)=j
     396  ! end if
     397  ! enddo
     398  ! do jj=1,num1
     399  ! j=index1(jj)
     400  ! write(6,*)  ' error = cloud fraction less than zero'
     401  ! write(6,*) ' or '
     402  ! write(6,*)  ' error = cloud fraction greater than 1'
     403  ! write(6,*) 'value at point ',j,' is ',cc(j,ilev)
     404  ! write(6,*) 'level ',ilev
     405  ! STOP
     406  ! enddo
     407  ! num1=0
     408  ! do j=1,npoints
     409  ! if (conv(j,ilev) .lt. 0. .or. conv(j,ilev) .gt. 1.) then
     410  ! num1=num1+1
     411  ! index1(num1)=j
     412  ! end if
     413  ! enddo
     414  ! do jj=1,num1
     415  ! j=index1(jj)
     416  ! write(6,*)
     417  ! &           ' error = convective cloud fraction less than zero'
     418  ! write(6,*) ' or '
     419  ! write(6,*)
     420  ! &           ' error = convective cloud fraction greater than 1'
     421  ! write(6,*) 'value at point ',j,' is ',conv(j,ilev)
     422  ! write(6,*) 'level ',ilev
     423  ! STOP
     424  ! enddo
     425
     426  ! num1=0
     427  ! do j=1,npoints
     428  ! if (dtau_s(j,ilev) .lt. 0.) then
     429  ! num1=num1+1
     430  ! index1(num1)=j
     431  ! end if
     432  ! enddo
     433  ! do jj=1,num1
     434  ! j=index1(jj)
     435  ! write(6,*)
     436  ! &           ' error = stratiform cloud opt. depth less than zero'
     437  ! write(6,*) 'value at point ',j,' is ',dtau_s(j,ilev)
     438  ! write(6,*) 'level ',ilev
     439  ! STOP
     440  ! enddo
     441  ! num1=0
     442  ! do j=1,npoints
     443  ! if (dtau_c(j,ilev) .lt. 0.) then
     444  ! num1=num1+1
     445  ! index1(num1)=j
     446  ! end if
     447  ! enddo
     448  ! do jj=1,num1
     449  ! j=index1(jj)
     450  ! write(6,*)
     451  ! &           ' error = convective cloud opt. depth less than zero'
     452  ! write(6,*) 'value at point ',j,' is ',dtau_c(j,ilev)
     453  ! write(6,*) 'level ',ilev
     454  ! STOP
     455  ! enddo
     456
     457  ! num1=0
     458  ! do j=1,npoints
     459  ! if (dem_s(j,ilev) .lt. 0. .or. dem_s(j,ilev) .gt. 1.) then
     460  ! num1=num1+1
     461  ! index1(num1)=j
     462  ! end if
     463  ! enddo
     464  ! do jj=1,num1
     465  ! j=index1(jj)
     466  ! write(6,*)
     467  ! &           ' error = stratiform cloud emissivity less than zero'
     468  ! write(6,*)'or'
     469  ! write(6,*)
     470  ! &           ' error = stratiform cloud emissivity greater than 1'
     471  ! write(6,*) 'value at point ',j,' is ',dem_s(j,ilev)
     472  ! write(6,*) 'level ',ilev
     473  ! STOP
     474  ! enddo
     475
     476  ! num1=0
     477  ! do j=1,npoints
     478  ! if (dem_c(j,ilev) .lt. 0. .or. dem_c(j,ilev) .gt. 1.) then
     479  ! num1=num1+1
     480  ! index1(num1)=j
     481  ! end if
     482  ! enddo
     483  ! do jj=1,num1
     484  ! j=index1(jj)
     485  ! write(6,*)
     486  ! &           ' error = convective cloud emissivity less than zero'
     487  ! write(6,*)'or'
     488  ! write(6,*)
     489  ! &           ' error = convective cloud emissivity greater than 1'
     490  ! write (6,*)
     491  ! &          'j=',j,'ilev=',ilev,'dem_c(j,ilev) =',dem_c(j,ilev)
     492  ! STOP
     493  ! enddo
     494  ! 13    continue
     495
     496
     497  DO ibox = 1, ncol
     498    DO j = 1, npoints
     499      boxpos(j, ibox) = (ibox-.5)/ncol
     500    END DO
     501  END DO
     502
     503  ! ---------------------------------------------------!
     504  ! Initialise working variables
     505  ! ---------------------------------------------------!
     506
     507  ! Initialised frac_out to zero
     508
     509  DO ilev = 1, nlev
     510    DO ibox = 1, ncol
     511      DO j = 1, npoints
     512        frac_out(j, ibox, ilev) = 0.0
     513      END DO
     514    END DO
     515  END DO
     516
     517  ! IM
     518  ! if (ncolprint.ne.0) then
     519  ! write (6,'(a)') 'frac_out_pp_rev:'
     520  ! do j=1,npoints,1000
     521  ! write(6,'(a10)') 'j='
     522  ! write(6,'(8I10)') j
     523  ! write (6,'(8f5.2)')
     524  ! &     ((frac_out(j,ibox,ilev),ibox=1,ncolprint),ilev=1,nlev)
     525
     526  ! enddo
     527  ! write (6,'(a)') 'ncol:'
     528  ! write (6,'(I3)') ncol
     529  ! endif
     530  ! if (ncolprint.ne.0) then
     531  ! write (6,'(a)') 'last_frac_pp:'
     532  ! do j=1,npoints,1000
     533  ! write(6,'(a10)') 'j='
     534  ! write(6,'(8I10)') j
     535  ! write (6,'(8f5.2)') (tca(j,0))
     536  ! enddo
     537  ! endif
     538
     539  ! ---------------------------------------------------!
     540  ! ALLOCATE CLOUD INTO BOXES, FOR NCOLUMNS, NLEVELS
     541  ! frac_out is the array that contains the information
     542  ! where 0 is no cloud, 1 is a stratiform cloud and 2 is a
     543  ! convective cloud
     544
     545    !loop over vertical levels
     546  DO ilev = 1, nlev
     547
     548    ! Initialise threshold
     549
     550    IF (ilev==1) THEN
     551        ! If max overlap
     552      IF (overlap==1) THEN
     553          ! select pixels spread evenly
     554          ! across the gridbox
     555        DO ibox = 1, ncol
     556          DO j = 1, npoints
     557            threshold(j, ibox) = boxpos(j, ibox)
     558          END DO
     559        END DO
     560      ELSE
     561        DO ibox = 1, ncol
     562          CALL ran0_vec(npoints, seed, ran)
     563            ! select random pixels from the non-convective
     564            ! part the gridbox ( some will be converted into
     565            ! convective pixels below )
     566          DO j = 1, npoints
     567            threshold(j, ibox) = cca(j, ilev) + (1-cca(j,ilev))*ran(j)
     568          END DO
     569        END DO
     570      END IF
     571      ! IM
     572      ! IF (ncolprint.ne.0) then
     573      ! write (6,'(a)') 'threshold_nsf2:'
     574      ! do j=1,npoints,1000
     575      ! write(6,'(a10)') 'j='
     576      ! write(6,'(8I10)') j
     577      ! write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint)
     578      ! enddo
     579      ! ENDIF
     580    END IF
     581
     582    ! IF (ncolprint.ne.0) then
     583    ! write (6,'(a)') 'ilev:'
     584    ! write (6,'(I2)') ilev
     585    ! ENDIF
     586
     587    DO ibox = 1, ncol
     588
     589        ! All versions
     590      DO j = 1, npoints
     591        IF (boxpos(j,ibox)<=cca(j,ilev)) THEN
     592          ! IM REAL           maxocc(j,ibox) = 1
     593          maxocc(j, ibox) = 1.0
     594        ELSE
     595          ! IM REAL           maxocc(j,ibox) = 0
     596          maxocc(j, ibox) = 0.0
     597        END IF
     598      END DO
     599
     600        ! Max overlap
     601      IF (overlap==1) THEN
     602        DO j = 1, npoints
     603          threshold_min(j, ibox) = cca(j, ilev)
     604          ! IM REAL           maxosc(j,ibox)=1
     605          maxosc(j, ibox) = 1.0
     606        END DO
     607      END IF
     608
     609        ! Random overlap
     610      IF (overlap==2) THEN
     611        DO j = 1, npoints
     612          threshold_min(j, ibox) = cca(j, ilev)
     613          ! IM REAL           maxosc(j,ibox)=0
     614          maxosc(j, ibox) = 0.0
     615        END DO
     616      END IF
     617
     618        ! Max/Random overlap
     619      IF (overlap==3) THEN
     620        DO j = 1, npoints
     621          threshold_min(j, ibox) = max(cca(j,ilev), min(tca(j,ilev-1),tca(j, &
     622            ilev)))
     623          IF (threshold(j,ibox)<min(tca(j,ilev-1),tca(j, &
     624              ilev)) .AND. (threshold(j,ibox)>cca(j,ilev))) THEN
     625            ! IM REAL                maxosc(j,ibox)= 1
     626            maxosc(j, ibox) = 1.0
     627          ELSE
     628            ! IM REAL                 maxosc(j,ibox)= 0
     629            maxosc(j, ibox) = 0.0
     630          END IF
     631        END DO
     632      END IF
     633
     634        ! Reset threshold
     635      CALL ran0_vec(npoints, seed, ran)
     636
     637      DO j = 1, npoints
     638        threshold(j, ibox) = &       !if max overlapped conv cloud
     639          maxocc(j, ibox)*(boxpos(j,ibox)) + &   !else
     640          (1-maxocc(j,ibox))*( &     !if max overlapped strat cloud
     641          (maxosc(j,ibox))*( &       !threshold=boxpos
     642          threshold(j,ibox))+ &      !else
     643          (1-maxosc(j,ibox))*( &     !threshold_min=random[thrmin,1]
     644          threshold_min(j,ibox)+(1-threshold_min(j,ibox))*ran(j)))
     645      END DO
     646
     647    END DO ! ibox
     648
     649    ! Fill frac_out with 1's where tca is greater than the threshold
     650
     651    DO ibox = 1, ncol
     652      DO j = 1, npoints
     653        IF (tca(j,ilev)>threshold(j,ibox)) THEN
     654          ! IM REAL             frac_out(j,ibox,ilev)=1
     655          frac_out(j, ibox, ilev) = 1.0
     656        ELSE
     657          ! IM REAL             frac_out(j,ibox,ilev)=0
     658          frac_out(j, ibox, ilev) = 0.0
     659        END IF
     660      END DO
     661    END DO
     662
     663    ! Code to partition boxes into startiform and convective parts
     664    ! goes here
     665
     666    DO ibox = 1, ncol
     667      DO j = 1, npoints
     668        IF (threshold(j,ibox)<=cca(j,ilev)) THEN
     669            ! = 2 IF threshold le cca(j)
     670          ! IM REAL                  frac_out(j,ibox,ilev) = 2
     671          frac_out(j, ibox, ilev) = 2.0
     672        ELSE
     673            ! = the same IF NOT threshold le cca(j)
     674          frac_out(j, ibox, ilev) = frac_out(j, ibox, ilev)
     675        END IF
     676      END DO
     677    END DO
     678
     679    ! Set last_frac to tca at this level, so as to be tca
     680    ! from last level next time round
     681
     682    ! IM
     683    ! if (ncolprint.ne.0) then
     684
     685    ! do j=1,npoints ,1000
     686    ! write(6,'(a10)') 'j='
     687    ! write(6,'(8I10)') j
     688    ! write (6,'(a)') 'last_frac:'
     689    ! write (6,'(8f5.2)') (tca(j,ilev-1))
     690
     691    ! write (6,'(a)') 'cca:'
     692    ! write (6,'(8f5.2)') (cca(j,ilev),ibox=1,ncolprint)
     693
     694    ! write (6,'(a)') 'max_overlap_cc:'
     695    ! write (6,'(8f5.2)') (maxocc(j,ibox),ibox=1,ncolprint)
     696
     697    ! write (6,'(a)') 'max_overlap_sc:'
     698    ! write (6,'(8f5.2)') (maxosc(j,ibox),ibox=1,ncolprint)
     699
     700    ! write (6,'(a)') 'threshold_min_nsf2:'
     701    ! write (6,'(8f5.2)') (threshold_min(j,ibox),ibox=1,ncolprint)
     702
     703    ! write (6,'(a)') 'threshold_nsf2:'
     704    ! write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint)
     705
     706    ! write (6,'(a)') 'frac_out_pp_rev:'
     707    ! write (6,'(8f5.2)')
     708    ! &       ((frac_out(j,ibox,ilev2),ibox=1,ncolprint),ilev2=1,nlev)
     709    ! enddo
     710    ! endif
     711
     712
     713  END DO
     714
     715  ! ---------------------------------------------------!
     716
     717
     718
     719  ! ---------------------------------------------------!
     720  ! COMPUTE CLOUD OPTICAL DEPTH FOR EACH COLUMN and
     721  ! put into vector tau
     722
     723    !initialize tau and albedocld to zero
     724  ! loop over nlev
     725  DO ibox = 1, ncol
     726    DO j = 1, npoints
     727      tau(j, ibox) = 0.
     728      albedocld(j, ibox) = 0.
     729      boxtau(j, ibox) = 0.
     730      boxptop(j, ibox) = 0.
     731      box_cloudy(j, ibox) = .FALSE.
     732    END DO
     733  END DO
     734
     735    !compute total cloud optical depth for each column
     736  DO ilev = 1, nlev
     737      !increment tau for each of the boxes
     738    DO ibox = 1, ncol
     739      DO j = 1, npoints
     740        ! IM REAL              if (frac_out(j,ibox,ilev).eq.1) then
     741        IF (frac_out(j,ibox,ilev)==1.0) THEN
     742          tau(j, ibox) = tau(j, ibox) + dtau_s(j, ilev)
     743        END IF
     744        ! IM REAL              if (frac_out(j,ibox,ilev).eq.2) then
     745        IF (frac_out(j,ibox,ilev)==2.0) THEN
     746          tau(j, ibox) = tau(j, ibox) + dtau_c(j, ilev)
     747        END IF
     748      END DO
     749    END DO ! ibox
     750  END DO ! ilev
     751  ! IM
     752  ! if (ncolprint.ne.0) then
     753
     754  ! do j=1,npoints ,1000
     755  ! write(6,'(a10)') 'j='
     756  ! write(6,'(8I10)') j
     757  ! write(6,'(i2,1X,8(f7.2,1X))')
     758  ! &          ilev,
     759  ! &          (tau(j,ibox),ibox=1,ncolprint)
     760  ! enddo
     761  ! endif
     762
     763  ! ---------------------------------------------------!
     764
     765
     766
     767
     768  ! ---------------------------------------------------!
     769  ! COMPUTE INFRARED BRIGHTNESS TEMPERUATRES
     770  ! AND CLOUD TOP TEMPERATURE SATELLITE SHOULD SEE
     771
     772  ! again this is only done if top_height = 1 or 3
     773
     774  ! fluxtop is the 10.5 micron radiance at the top of the
     775  ! atmosphere
     776  ! trans_layers_above is the total transmissivity in the layers
     777  ! above the current layer
     778  ! fluxtop_clrsky(j) and trans_layers_above_clrsky(j) are the clear
     779  ! sky versions of these quantities.
     780
     781  IF (top_height==1 .OR. top_height==3) THEN
     782
     783
     784      !----------------------------------------------------------------------
     785      !
     786      !             DO CLEAR SKY RADIANCE CALCULATION FIRST
     787      !
     788      !compute water vapor continuum emissivity
     789      !this treatment follows Schwarkzopf and Ramasamy
     790      !JGR 1999,vol 104, pages 9467-9499.
     791      !the emissivity is calculated at a wavenumber of 955 cm-1,
     792      !or 10.47 microns
     793    wtmair = 28.9644
     794    wtmh20 = 18.01534
     795    navo = 6.023E+23
     796    grav = 9.806650E+02
     797    pstd = 1.013250E+06
     798    t0 = 296.
     799    ! IM
     800    ! if (ncolprint .ne. 0)
     801    ! &         write(6,*)  'ilev   pw (kg/m2)   tauwv(j)      dem_wv'
     802    DO ilev = 1, nlev
     803      DO j = 1, npoints
     804          !press and dpress are dyne/cm2 = Pascals *10
     805        press(j) = pfull(j, ilev)*10.
     806        dpress(j) = (phalf(j,ilev+1)-phalf(j,ilev))*10
     807          !atmden = g/cm2 = kg/m2 / 10
     808        atmden(j) = dpress(j)/grav
     809        rvh20(j) = qv(j, ilev)*wtmair/wtmh20
     810        wk(j) = rvh20(j)*navo*atmden(j)/wtmair
     811        rhoave(j) = (press(j)/pstd)*(t0/at(j,ilev))
     812        rh20s(j) = rvh20(j)*rhoave(j)
     813        rfrgn(j) = rhoave(j) - rh20s(j)
     814        tmpexp(j) = exp(-0.02*(at(j,ilev)-t0))
     815        tauwv(j) = wk(j)*1.E-20*((0.0224697*rh20s(j)*tmpexp(j))+(3.41817E-7* &
     816          rfrgn(j)))*0.98
     817        dem_wv(j, ilev) = 1. - exp(-1.*tauwv(j))
     818      END DO
     819      ! IM
     820      ! if (ncolprint .ne. 0) then
     821      ! do j=1,npoints ,1000
     822      ! write(6,'(a10)') 'j='
     823      ! write(6,'(8I10)') j
     824      ! write(6,'(i2,1X,3(f8.3,3X))') ilev,
     825      ! &           qv(j,ilev)*(phalf(j,ilev+1)-phalf(j,ilev))/(grav/100.),
     826      ! &           tauwv(j),dem_wv(j,ilev)
     827      ! enddo
     828      ! endif
     829    END DO
     830
     831      !initialize variables
     832    DO j = 1, npoints
     833      fluxtop_clrsky(j) = 0.
     834      trans_layers_above_clrsky(j) = 1.
     835    END DO
     836
     837    DO ilev = 1, nlev
     838      DO j = 1, npoints
     839
     840          ! Black body emission at temperature of the layer
     841
     842        bb(j) = 1/(exp(1307.27/at(j,ilev))-1.)
     843          !bb(j)= 5.67e-8*at(j,ilev)**4
     844
     845          ! increase TOA flux by flux emitted from layer
     846          ! times total transmittance in layers above
     847
     848        fluxtop_clrsky(j) = fluxtop_clrsky(j) + dem_wv(j, ilev)*bb(j)* &
     849          trans_layers_above_clrsky(j)
     850
     851          ! update trans_layers_above with transmissivity
     852          ! from this layer for next time around loop
     853
     854        trans_layers_above_clrsky(j) = trans_layers_above_clrsky(j)* &
     855          (1.-dem_wv(j,ilev))
     856
     857
     858      END DO
     859      ! IM
     860      ! if (ncolprint.ne.0) then
     861      ! do j=1,npoints ,1000
     862      ! write(6,'(a10)') 'j='
     863      ! write(6,'(8I10)') j
     864      ! write (6,'(a)') 'ilev:'
     865      ! write (6,'(I2)') ilev
     866
     867      ! write (6,'(a)')
     868      ! &        'emiss_layer,100.*bb(j),100.*f,total_trans:'
     869      ! write (6,'(4(f7.2,1X))') dem_wv(j,ilev),100.*bb(j),
     870      ! &             100.*fluxtop_clrsky(j),trans_layers_above_clrsky(j)
     871      ! enddo
     872      ! endif
     873
     874    END DO !loop over level
     875
     876    DO j = 1, npoints
     877        !add in surface emission
     878      bb(j) = 1/(exp(1307.27/skt(j))-1.)
     879        !bb(j)=5.67e-8*skt(j)**4
     880
     881      fluxtop_clrsky(j) = fluxtop_clrsky(j) + emsfc_lw*bb(j)* &
     882        trans_layers_above_clrsky(j)
     883    END DO
     884
     885    ! IM
     886    ! if (ncolprint.ne.0) then
     887    ! do j=1,npoints ,1000
     888    ! write(6,'(a10)') 'j='
     889    ! write(6,'(8I10)') j
     890    ! write (6,'(a)') 'id:'
     891    ! write (6,'(a)') 'surface'
     892
     893    ! write (6,'(a)') 'emsfc,100.*bb(j),100.*f,total_trans:'
     894    ! write (6,'(4(f7.2,1X))') emsfc_lw,100.*bb(j),
     895    ! &      100.*fluxtop_clrsky(j),
     896    ! &       trans_layers_above_clrsky(j)
     897    ! enddo
     898    ! endif
     899
     900
     901      !
     902      !           END OF CLEAR SKY CALCULATION
     903      !
     904      !----------------------------------------------------------------
     905
     906
     907    ! IM
     908    ! if (ncolprint.ne.0) then
     909
     910    ! do j=1,npoints ,1000
     911    ! write(6,'(a10)') 'j='
     912    ! write(6,'(8I10)') j
     913    ! write (6,'(a)') 'ts:'
     914    ! write (6,'(8f7.2)') (skt(j),ibox=1,ncolprint)
     915
     916    ! write (6,'(a)') 'ta_rev:'
     917    ! write (6,'(8f7.2)')
     918    ! &       ((at(j,ilev2),ibox=1,ncolprint),ilev2=1,nlev)
     919
     920    ! enddo
     921    ! endif
     922      !loop over columns
     923    DO ibox = 1, ncol
     924      DO j = 1, npoints
     925        fluxtop(j, ibox) = 0.
     926        trans_layers_above(j, ibox) = 1.
     927      END DO
     928    END DO
     929
     930    DO ilev = 1, nlev
     931      DO j = 1, npoints
     932          ! Black body emission at temperature of the layer
     933
     934        bb(j) = 1/(exp(1307.27/at(j,ilev))-1.)
     935          !bb(j)= 5.67e-8*at(j,ilev)**4
     936      END DO
     937
     938      DO ibox = 1, ncol
     939        DO j = 1, npoints
     940
     941            ! emissivity for point in this layer
     942          ! IM REAL             if (frac_out(j,ibox,ilev).eq.1) then
     943          IF (frac_out(j,ibox,ilev)==1.0) THEN
     944            dem(j, ibox) = 1. - ((1.-dem_wv(j,ilev))*(1.-dem_s(j,ilev)))
     945            ! IM REAL             else if (frac_out(j,ibox,ilev).eq.2) then
     946          ELSE IF (frac_out(j,ibox,ilev)==2.0) THEN
     947            dem(j, ibox) = 1. - ((1.-dem_wv(j,ilev))*(1.-dem_c(j,ilev)))
     948          ELSE
     949            dem(j, ibox) = dem_wv(j, ilev)
     950          END IF
     951
     952
     953            ! increase TOA flux by flux emitted from layer
     954            ! times total transmittance in layers above
     955
     956          fluxtop(j, ibox) = fluxtop(j, ibox) + dem(j, ibox)*bb(j)* &
     957            trans_layers_above(j, ibox)
     958
     959            ! update trans_layers_above with transmissivity
     960            ! from this layer for next time around loop
     961
     962          trans_layers_above(j, ibox) = trans_layers_above(j, ibox)* &
     963            (1.-dem(j,ibox))
     964
     965        END DO ! j
     966      END DO ! ibox
     967
     968      ! IM
     969      ! if (ncolprint.ne.0) then
     970      ! do j=1,npoints,1000
     971      ! write (6,'(a)') 'ilev:'
     972      ! write (6,'(I2)') ilev
     973
     974      ! write(6,'(a10)') 'j='
     975      ! write(6,'(8I10)') j
     976      ! write (6,'(a)') 'emiss_layer:'
     977      ! write (6,'(8f7.2)') (dem(j,ibox),ibox=1,ncolprint)
     978
     979      ! write (6,'(a)') '100.*bb(j):'
     980      ! write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint)
     981
     982      ! write (6,'(a)') '100.*f:'
     983      ! write (6,'(8f7.2)')
     984      ! &         (100.*fluxtop(j,ibox),ibox=1,ncolprint)
     985
     986      ! write (6,'(a)') 'total_trans:'
     987      ! write (6,'(8f7.2)')
     988      ! &          (trans_layers_above(j,ibox),ibox=1,ncolprint)
     989      ! enddo
     990      ! endif
     991
     992    END DO ! ilev
     993
     994
     995    DO j = 1, npoints
     996        !add in surface emission
     997      bb(j) = 1/(exp(1307.27/skt(j))-1.)
     998        !bb(j)=5.67e-8*skt(j)**4
     999    END DO
     1000
     1001    DO ibox = 1, ncol
     1002      DO j = 1, npoints
     1003
     1004          !add in surface emission
     1005
     1006        fluxtop(j, ibox) = fluxtop(j, ibox) + emsfc_lw*bb(j)* &
     1007          trans_layers_above(j, ibox)
     1008
     1009      END DO
     1010    END DO
     1011
     1012    ! IM
     1013    ! if (ncolprint.ne.0) then
     1014
     1015    ! do j=1,npoints ,1000
     1016    ! write(6,'(a10)') 'j='
     1017    ! write(6,'(8I10)') j
     1018    ! write (6,'(a)') 'id:'
     1019    ! write (6,'(a)') 'surface'
     1020
     1021    ! write (6,'(a)') 'emiss_layer:'
     1022    ! write (6,'(8f7.2)') (dem(1,ibox),ibox=1,ncolprint)
     1023
     1024    ! write (6,'(a)') '100.*bb(j):'
     1025    ! write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint)
     1026
     1027    ! write (6,'(a)') '100.*f:'
     1028    ! write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint)
     1029    ! end do
     1030    ! endif
     1031
     1032      !now that you have the top of atmosphere radiance account
     1033      !for ISCCP procedures to determine cloud top temperature
     1034
     1035      !account for partially transmitting cloud recompute flux
     1036      !ISCCP would see assuming a single layer cloud
     1037      !note choice here of 2.13, as it is primarily ice
     1038      !clouds which have partial emissivity and need the
     1039      !adjustment performed in this section
     1040      !
     1041      !If it turns out that the cloud brightness temperature
     1042      !is greater than 260K, then the liquid cloud conversion
     1043      !factor of 2.56 is used.
     1044      !
     1045      !Note that this is discussed on pages 85-87 of
     1046      !the ISCCP D level documentation (Rossow et al. 1996)
     1047
     1048    DO j = 1, npoints
     1049        !compute minimum brightness temperature and optical depth
     1050      btcmin(j) = 1./(exp(1307.27/(attrop(j)-5.))-1.)
     1051    END DO
     1052    DO ibox = 1, ncol
     1053      DO j = 1, npoints
     1054        transmax(j) = (fluxtop(j,ibox)-btcmin(j))/(fluxtop_clrsky(j)-btcmin(j &
     1055          ))
     1056          !note that the initial setting of tauir(j) is needed so that
     1057          !tauir(j) has a realistic value should the next if block be
     1058          !bypassed
     1059        tauir(j) = tau(j, ibox)*rec2p13
     1060        taumin(j) = -1.*log(max(min(transmax(j),0.9999999),0.001))
     1061
     1062      END DO
     1063
     1064      IF (top_height==1) THEN
     1065        DO j = 1, npoints
     1066          IF (transmax(j)>0.001 .AND. transmax(j)<=0.9999999) THEN
     1067            fluxtopinit(j) = fluxtop(j, ibox)
     1068            tauir(j) = tau(j, ibox)*rec2p13
     1069          END IF
     1070        END DO
     1071        DO icycle = 1, 2
     1072          DO j = 1, npoints
     1073            IF (tau(j,ibox)>(tauchk)) THEN
     1074              IF (transmax(j)>0.001 .AND. transmax(j)<=0.9999999) THEN
     1075                emcld(j, ibox) = 1. - exp(-1.*tauir(j))
     1076                fluxtop(j, ibox) = fluxtopinit(j) - ((1.-emcld(j, &
     1077                  ibox))*fluxtop_clrsky(j))
     1078                fluxtop(j, ibox) = max(1.E-06, (fluxtop(j,ibox)/emcld(j, &
     1079                  ibox)))
     1080                tb(j, ibox) = 1307.27/(log(1.+(1./fluxtop(j,ibox))))
     1081                IF (tb(j,ibox)>260.) THEN
     1082                  tauir(j) = tau(j, ibox)/2.56
     1083                END IF
     1084              END IF
     1085            END IF
     1086          END DO
     1087        END DO
     1088
     1089      END IF
     1090
     1091      DO j = 1, npoints
     1092        IF (tau(j,ibox)>(tauchk)) THEN
     1093            !cloudy box
     1094          tb(j, ibox) = 1307.27/(log(1.+(1./fluxtop(j,ibox))))
     1095          IF (top_height==1 .AND. tauir(j)<taumin(j)) THEN
     1096            tb(j, ibox) = attrop(j) - 5.
     1097            tau(j, ibox) = 2.13*taumin(j)
     1098          END IF
     1099        ELSE
     1100            !clear sky brightness temperature
     1101          tb(j, ibox) = 1307.27/(log(1.+(1./fluxtop_clrsky(j))))
     1102        END IF
     1103      END DO ! j
     1104    END DO ! ibox
     1105
     1106    ! IM
     1107    ! if (ncolprint.ne.0) then
     1108
     1109    ! do j=1,npoints,1000
     1110    ! write(6,'(a10)') 'j='
     1111    ! write(6,'(8I10)') j
     1112
     1113    ! write (6,'(a)') 'attrop:'
     1114    ! write (6,'(8f7.2)') (attrop(j))
     1115
     1116    ! write (6,'(a)') 'btcmin:'
     1117    ! write (6,'(8f7.2)') (btcmin(j))
     1118
     1119    ! write (6,'(a)') 'fluxtop_clrsky*100:'
     1120    ! write (6,'(8f7.2)')
     1121    ! &      (100.*fluxtop_clrsky(j))
     1122
     1123    ! write (6,'(a)') '100.*f_adj:'
     1124    ! write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint)
     1125
     1126    ! write (6,'(a)') 'transmax:'
     1127    ! write (6,'(8f7.2)') (transmax(ibox),ibox=1,ncolprint)
     1128
     1129    ! write (6,'(a)') 'tau:'
     1130    ! write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint)
     1131
     1132    ! write (6,'(a)') 'emcld:'
     1133    ! write (6,'(8f7.2)') (emcld(j,ibox),ibox=1,ncolprint)
     1134
     1135    ! write (6,'(a)') 'total_trans:'
     1136    ! write (6,'(8f7.2)')
     1137    ! &          (trans_layers_above(j,ibox),ibox=1,ncolprint)
     1138
     1139    ! write (6,'(a)') 'total_emiss:'
     1140    ! write (6,'(8f7.2)')
     1141    ! &          (1.0-trans_layers_above(j,ibox),ibox=1,ncolprint)
     1142
     1143    ! write (6,'(a)') 'total_trans:'
     1144    ! write (6,'(8f7.2)')
     1145    ! &          (trans_layers_above(j,ibox),ibox=1,ncolprint)
     1146
     1147    ! write (6,'(a)') 'ppout:'
     1148    ! write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint)
     1149    ! enddo ! j
     1150    ! endif
     1151
     1152  END IF
     1153
     1154  ! ---------------------------------------------------!
     1155
     1156
     1157  ! ---------------------------------------------------!
     1158  ! DETERMINE CLOUD TOP PRESSURE
     1159
     1160  ! again the 2 methods differ according to whether
     1161  ! or not you use the physical cloud top pressure (top_height = 2)
     1162  ! or the radiatively determined cloud top pressure (top_height = 1 or 3)
     1163
     1164
     1165    !compute cloud top pressure
     1166  DO ibox = 1, ncol
     1167      !segregate according to optical thickness
     1168    IF (top_height==1 .OR. top_height==3) THEN
     1169        !find level whose temperature
     1170        !most closely matches brightness temperature
     1171      DO j = 1, npoints
     1172        nmatch(j) = 0
     1173      END DO
     1174      DO ilev = 1, nlev - 1
     1175        ! cdir nodep
     1176        DO j = 1, npoints
     1177          IF ((at(j,ilev)>=tb(j,ibox) .AND. at(j,ilev+1)<tb(j, &
     1178              ibox)) .OR. (at(j,ilev)<=tb(j,ibox) .AND. at(j,ilev+1)>tb(j, &
     1179              ibox))) THEN
     1180
     1181            nmatch(j) = nmatch(j) + 1
     1182            IF (abs(at(j,ilev)-tb(j,ibox))<abs(at(j,ilev+1)-tb(j,ibox))) THEN
     1183              match(j, nmatch(j)) = ilev
    5831184            ELSE
    584               DO ibox=1,ncol
    585                 call ran0_vec(npoints,seed,ran)
    586                 ! select random pixels from the non-convective
    587                 ! part the gridbox ( some will be converted into
    588                 ! convective pixels below )
    589                 do j=1,npoints
    590                   threshold(j,ibox)=
    591      &            cca(j,ilev)+(1-cca(j,ilev))*ran(j)
    592                 enddo
    593               enddo
    594             ENDIF
    595 cIM
    596 c            IF (ncolprint.ne.0) then
    597 c              write (6,'(a)') 'threshold_nsf2:'
    598 c                do j=1,npoints,1000
    599 c                write(6,'(a10)') 'j='
    600 c                write(6,'(8I10)') j
    601 c                write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint)
    602 c                enddo
    603 c            ENDIF
    604         ENDIF
    605 
    606 c        IF (ncolprint.ne.0) then
    607 c            write (6,'(a)') 'ilev:'
    608 c            write (6,'(I2)') ilev
    609 c        ENDIF
    610 
    611         DO ibox=1,ncol
    612 
    613           ! All versions
    614           do j=1,npoints
    615             if (boxpos(j,ibox).le.cca(j,ilev)) then
    616 cIM REAL           maxocc(j,ibox) = 1
    617               maxocc(j,ibox) = 1.0
    618             else
    619 cIM REAL           maxocc(j,ibox) = 0
    620               maxocc(j,ibox) = 0.0
    621             end if
    622           enddo
    623 
    624           ! Max overlap
    625           if (overlap.eq.1) then
    626             do j=1,npoints
    627               threshold_min(j,ibox)=cca(j,ilev)
    628 cIM REAL           maxosc(j,ibox)=1
    629               maxosc(j,ibox)=1.0
    630             enddo
    631           endif
    632 
    633           ! Random overlap
    634           if (overlap.eq.2) then
    635             do j=1,npoints
    636               threshold_min(j,ibox)=cca(j,ilev)
    637 cIM REAL           maxosc(j,ibox)=0
    638               maxosc(j,ibox)=0.0
    639             enddo
    640           endif
    641 
    642           ! Max/Random overlap
    643           if (overlap.eq.3) then
    644             do j=1,npoints
    645               threshold_min(j,ibox)=max(cca(j,ilev),
    646      &          min(tca(j,ilev-1),tca(j,ilev)))
    647               if (threshold(j,ibox)
    648      &          .lt.min(tca(j,ilev-1),tca(j,ilev))
    649      &          .and.(threshold(j,ibox).gt.cca(j,ilev))) then
    650 cIM REAL                maxosc(j,ibox)= 1
    651                    maxosc(j,ibox)= 1.0
    652               else
    653 cIM REAL                 maxosc(j,ibox)= 0
    654                    maxosc(j,ibox)= 0.0
    655               end if
    656             enddo
    657           endif
    658    
    659           ! Reset threshold
    660           call ran0_vec(npoints,seed,ran)
    661            
    662           do j=1,npoints
    663             threshold(j,ibox)=
    664               !if max overlapped conv cloud
    665      &        maxocc(j,ibox) * (                                       
    666      &            boxpos(j,ibox)                                               
    667      &        ) +                                                     
    668               !else
    669      &        (1-maxocc(j,ibox)) * (                                   
    670                   !if max overlapped strat cloud
    671      &            (maxosc(j,ibox)) * (                                 
    672                       !threshold=boxpos
    673      &                threshold(j,ibox)                                       
    674      &            ) +                                                 
    675                   !else
    676      &            (1-maxosc(j,ibox)) * (                               
    677                       !threshold_min=random[thrmin,1]
    678      &                threshold_min(j,ibox)+
    679      &                  (1-threshold_min(j,ibox))*ran(j) 
    680      &           )
    681      &        )
    682           enddo
    683 
    684         ENDDO ! ibox
    685 
    686 !          Fill frac_out with 1's where tca is greater than the threshold
    687 
    688            DO ibox=1,ncol
    689              do j=1,npoints
    690                if (tca(j,ilev).gt.threshold(j,ibox)) then
    691 cIM REAL             frac_out(j,ibox,ilev)=1
    692                frac_out(j,ibox,ilev)=1.0
    693                else
    694 cIM REAL             frac_out(j,ibox,ilev)=0
    695                frac_out(j,ibox,ilev)=0.0
    696                end if               
    697              enddo
    698            ENDDO
    699 
    700 !           Code to partition boxes into startiform and convective parts
    701 !           goes here
    702 
    703            DO ibox=1,ncol
    704              do j=1,npoints
    705                 if (threshold(j,ibox).le.cca(j,ilev)) then
    706                     ! = 2 IF threshold le cca(j)
    707 cIM REAL                  frac_out(j,ibox,ilev) = 2
    708                     frac_out(j,ibox,ilev) = 2.0
    709                 else
    710                     ! = the same IF NOT threshold le cca(j)
    711                     frac_out(j,ibox,ilev) = frac_out(j,ibox,ilev)
    712                 end if
    713              enddo
    714            ENDDO
    715 
    716 !         Set last_frac to tca at this level, so as to be tca
    717 !         from last level next time round
    718 
    719 cIM
    720 c          if (ncolprint.ne.0) then
    721 
    722 c            do j=1,npoints ,1000
    723 c            write(6,'(a10)') 'j='
    724 c            write(6,'(8I10)') j
    725 c            write (6,'(a)') 'last_frac:'
    726 c            write (6,'(8f5.2)') (tca(j,ilev-1))
    727    
    728 c            write (6,'(a)') 'cca:'
    729 c            write (6,'(8f5.2)') (cca(j,ilev),ibox=1,ncolprint)
    730    
    731 c            write (6,'(a)') 'max_overlap_cc:'
    732 c            write (6,'(8f5.2)') (maxocc(j,ibox),ibox=1,ncolprint)
    733    
    734 c            write (6,'(a)') 'max_overlap_sc:'
    735 c            write (6,'(8f5.2)') (maxosc(j,ibox),ibox=1,ncolprint)
    736    
    737 c            write (6,'(a)') 'threshold_min_nsf2:'
    738 c            write (6,'(8f5.2)') (threshold_min(j,ibox),ibox=1,ncolprint)
    739    
    740 c            write (6,'(a)') 'threshold_nsf2:'
    741 c            write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint)
    742    
    743 c            write (6,'(a)') 'frac_out_pp_rev:'
    744 c            write (6,'(8f5.2)')
    745 c     &       ((frac_out(j,ibox,ilev2),ibox=1,ncolprint),ilev2=1,nlev)
    746 c            enddo
    747 c          endif
    748 
    749 200   CONTINUE    !loop over nlev
    750 
    751 !
    752 !     ---------------------------------------------------!
    753 
    754      
    755 !
    756 !     ---------------------------------------------------!
    757 !     COMPUTE CLOUD OPTICAL DEPTH FOR EACH COLUMN and
    758 !     put into vector tau
    759  
    760       !initialize tau and albedocld to zero
    761       do 15 ibox=1,ncol
    762         do j=1,npoints
    763             tau(j,ibox)=0.
    764             albedocld(j,ibox)=0.
    765             boxtau(j,ibox)=0.
    766             boxptop(j,ibox)=0.
    767             box_cloudy(j,ibox)=.false.
    768         enddo
    769 15    continue
    770 
    771       !compute total cloud optical depth for each column     
    772       do ilev=1,nlev
    773             !increment tau for each of the boxes
    774             do ibox=1,ncol
    775               do j=1,npoints
    776 cIM REAL              if (frac_out(j,ibox,ilev).eq.1) then
    777                  if (frac_out(j,ibox,ilev).eq.1.0) then
    778                         tau(j,ibox)=tau(j,ibox)
    779      &                     + dtau_s(j,ilev)
    780                  endif
    781 cIM REAL              if (frac_out(j,ibox,ilev).eq.2) then
    782                  if (frac_out(j,ibox,ilev).eq.2.0) then
    783                         tau(j,ibox)=tau(j,ibox)
    784      &                     + dtau_c(j,ilev)
    785                  end if
    786               enddo
    787             enddo ! ibox
    788       enddo ! ilev
    789 cIM
    790 c          if (ncolprint.ne.0) then
    791 
    792 c              do j=1,npoints ,1000
    793 c                write(6,'(a10)') 'j='
    794 c                write(6,'(8I10)') j
    795 c                write(6,'(i2,1X,8(f7.2,1X))')
    796 c     &          ilev,
    797 c     &          (tau(j,ibox),ibox=1,ncolprint)
    798 c              enddo
    799 c          endif
    800 !
    801 !     ---------------------------------------------------!
    802 
    803 
    804 
    805 !     
    806 !     ---------------------------------------------------!
    807 !     COMPUTE INFRARED BRIGHTNESS TEMPERUATRES
    808 !     AND CLOUD TOP TEMPERATURE SATELLITE SHOULD SEE
    809 !
    810 !     again this is only done if top_height = 1 or 3
    811 !
    812 !     fluxtop is the 10.5 micron radiance at the top of the
    813 !              atmosphere
    814 !     trans_layers_above is the total transmissivity in the layers
    815 !             above the current layer
    816 !     fluxtop_clrsky(j) and trans_layers_above_clrsky(j) are the clear
    817 !             sky versions of these quantities.
    818 
    819       if (top_height .eq. 1 .or. top_height .eq. 3) then
    820 
    821 
    822         !----------------------------------------------------------------------
    823         !   
    824         !             DO CLEAR SKY RADIANCE CALCULATION FIRST
    825         !
    826         !compute water vapor continuum emissivity
    827         !this treatment follows Schwarkzopf and Ramasamy
    828         !JGR 1999,vol 104, pages 9467-9499.
    829         !the emissivity is calculated at a wavenumber of 955 cm-1,
    830         !or 10.47 microns
    831         wtmair = 28.9644
    832         wtmh20 = 18.01534
    833         Navo = 6.023E+23
    834         grav = 9.806650E+02
    835         pstd = 1.013250E+06
    836         t0 = 296.
    837 cIM
    838 c        if (ncolprint .ne. 0)
    839 c     &         write(6,*)  'ilev   pw (kg/m2)   tauwv(j)      dem_wv'
    840         do 125 ilev=1,nlev
    841           do j=1,npoints
    842                !press and dpress are dyne/cm2 = Pascals *10
    843                press(j) = pfull(j,ilev)*10.
    844                dpress(j) = (phalf(j,ilev+1)-phalf(j,ilev))*10
    845                !atmden = g/cm2 = kg/m2 / 10
    846                atmden(j) = dpress(j)/grav
    847                rvh20(j) = qv(j,ilev)*wtmair/wtmh20
    848                wk(j) = rvh20(j)*Navo*atmden(j)/wtmair
    849                rhoave(j) = (press(j)/pstd)*(t0/at(j,ilev))
    850                rh20s(j) = rvh20(j)*rhoave(j)
    851                rfrgn(j) = rhoave(j)-rh20s(j)
    852                tmpexp(j) = exp(-0.02*(at(j,ilev)-t0))
    853                tauwv(j) = wk(j)*1.e-20*(
    854      &           (0.0224697*rh20s(j)*tmpexp(j)) +
    855      &                (3.41817e-7*rfrgn(j)) )*0.98
    856                dem_wv(j,ilev) = 1. - exp( -1. * tauwv(j))
    857           enddo
    858 cIM
    859 c               if (ncolprint .ne. 0) then
    860 c               do j=1,npoints ,1000
    861 c               write(6,'(a10)') 'j='
    862 c               write(6,'(8I10)') j
    863 c               write(6,'(i2,1X,3(f8.3,3X))') ilev,
    864 c     &           qv(j,ilev)*(phalf(j,ilev+1)-phalf(j,ilev))/(grav/100.),
    865 c     &           tauwv(j),dem_wv(j,ilev)
    866 c               enddo
    867 c               endif
    868 125     continue
    869 
    870         !initialize variables
    871         do j=1,npoints
    872           fluxtop_clrsky(j) = 0.
    873           trans_layers_above_clrsky(j)=1.
    874         enddo
    875 
    876         do ilev=1,nlev
    877           do j=1,npoints
    878  
    879             ! Black body emission at temperature of the layer
    880 
    881                 bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. )
    882                 !bb(j)= 5.67e-8*at(j,ilev)**4
    883 
    884                 ! increase TOA flux by flux emitted from layer
    885                 ! times total transmittance in layers above
    886 
    887                 fluxtop_clrsky(j) = fluxtop_clrsky(j)
    888      &            + dem_wv(j,ilev)*bb(j)*trans_layers_above_clrsky(j)
    889            
    890                 ! update trans_layers_above with transmissivity
    891                 ! from this layer for next time around loop
    892 
    893                 trans_layers_above_clrsky(j)=
    894      &            trans_layers_above_clrsky(j)*(1.-dem_wv(j,ilev))
    895                    
    896 
    897           enddo
    898 cIM 
    899 c            if (ncolprint.ne.0) then
    900 c             do j=1,npoints ,1000
    901 c              write(6,'(a10)') 'j='
    902 c              write(6,'(8I10)') j
    903 c              write (6,'(a)') 'ilev:'
    904 c              write (6,'(I2)') ilev
    905    
    906 c              write (6,'(a)')
    907 c     &        'emiss_layer,100.*bb(j),100.*f,total_trans:'
    908 c              write (6,'(4(f7.2,1X))') dem_wv(j,ilev),100.*bb(j),
    909 c     &             100.*fluxtop_clrsky(j),trans_layers_above_clrsky(j)
    910 c             enddo   
    911 c            endif
    912 
    913         enddo   !loop over level
    914        
    915         do j=1,npoints
    916           !add in surface emission
    917           bb(j)=1/( exp(1307.27/skt(j)) - 1. )
    918           !bb(j)=5.67e-8*skt(j)**4
    919 
    920           fluxtop_clrsky(j) = fluxtop_clrsky(j) + emsfc_lw * bb(j)
    921      &     * trans_layers_above_clrsky(j)
    922         enddo
    923 
    924 cIM
    925 c        if (ncolprint.ne.0) then
    926 c        do j=1,npoints ,1000
    927 c          write(6,'(a10)') 'j='
    928 c          write(6,'(8I10)') j
    929 c          write (6,'(a)') 'id:'
    930 c          write (6,'(a)') 'surface'
    931 
    932 c          write (6,'(a)') 'emsfc,100.*bb(j),100.*f,total_trans:'
    933 c          write (6,'(4(f7.2,1X))') emsfc_lw,100.*bb(j),
    934 c     &      100.*fluxtop_clrsky(j),
    935 c     &       trans_layers_above_clrsky(j)
    936 c        enddo
    937 c        endif
    938    
    939 
    940         !
    941         !           END OF CLEAR SKY CALCULATION
    942         !
    943         !----------------------------------------------------------------
    944 
    945 
    946 cIM
    947 c        if (ncolprint.ne.0) then
    948 
    949 c        do j=1,npoints ,1000
    950 c            write(6,'(a10)') 'j='
    951 c            write(6,'(8I10)') j
    952 c            write (6,'(a)') 'ts:'
    953 c            write (6,'(8f7.2)') (skt(j),ibox=1,ncolprint)
    954    
    955 c            write (6,'(a)') 'ta_rev:'
    956 c            write (6,'(8f7.2)')
    957 c     &       ((at(j,ilev2),ibox=1,ncolprint),ilev2=1,nlev)
    958 
    959 c        enddo
    960 c        endif
    961         !loop over columns
    962         do ibox=1,ncol
    963           do j=1,npoints
    964             fluxtop(j,ibox)=0.
    965             trans_layers_above(j,ibox)=1.
    966           enddo
    967         enddo
    968 
    969         do ilev=1,nlev
    970               do j=1,npoints
    971                 ! Black body emission at temperature of the layer
    972 
    973                 bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. )
    974                 !bb(j)= 5.67e-8*at(j,ilev)**4
    975               enddo
    976 
    977             do ibox=1,ncol
    978               do j=1,npoints
    979 
    980                 ! emissivity for point in this layer
    981 cIM REAL             if (frac_out(j,ibox,ilev).eq.1) then
    982                 if (frac_out(j,ibox,ilev).eq.1.0) then
    983                 dem(j,ibox)= 1. -
    984      &             ( (1. - dem_wv(j,ilev)) * (1. -  dem_s(j,ilev)) )
    985 cIM REAL             else if (frac_out(j,ibox,ilev).eq.2) then
    986                 else if (frac_out(j,ibox,ilev).eq.2.0) then
    987                 dem(j,ibox)= 1. -
    988      &             ( (1. - dem_wv(j,ilev)) * (1. -  dem_c(j,ilev)) )
    989                 else
    990                 dem(j,ibox)=  dem_wv(j,ilev)
    991                 end if
    992                
    993 
    994                 ! increase TOA flux by flux emitted from layer
    995                 ! times total transmittance in layers above
    996 
    997                 fluxtop(j,ibox) = fluxtop(j,ibox)
    998      &            + dem(j,ibox) * bb(j)
    999      &            * trans_layers_above(j,ibox)
    1000            
    1001                 ! update trans_layers_above with transmissivity
    1002                 ! from this layer for next time around loop
    1003 
    1004                 trans_layers_above(j,ibox)=
    1005      &            trans_layers_above(j,ibox)*(1.-dem(j,ibox))
    1006 
    1007               enddo ! j
    1008             enddo ! ibox
    1009 
    1010 cIM
    1011 c            if (ncolprint.ne.0) then
    1012 c              do j=1,npoints,1000
    1013 c              write (6,'(a)') 'ilev:'
    1014 c              write (6,'(I2)') ilev
    1015    
    1016 c              write(6,'(a10)') 'j='
    1017 c              write(6,'(8I10)') j
    1018 c              write (6,'(a)') 'emiss_layer:'
    1019 c              write (6,'(8f7.2)') (dem(j,ibox),ibox=1,ncolprint)
    1020        
    1021 c              write (6,'(a)') '100.*bb(j):'
    1022 c              write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint)
    1023        
    1024 c              write (6,'(a)') '100.*f:'
    1025 c              write (6,'(8f7.2)')
    1026 c     &         (100.*fluxtop(j,ibox),ibox=1,ncolprint)
    1027        
    1028 c              write (6,'(a)') 'total_trans:'
    1029 c              write (6,'(8f7.2)')
    1030 c     &          (trans_layers_above(j,ibox),ibox=1,ncolprint)
    1031 c              enddo
    1032 c          endif
    1033 
    1034         enddo ! ilev
    1035 
    1036 
    1037           do j=1,npoints
    1038             !add in surface emission
    1039             bb(j)=1/( exp(1307.27/skt(j)) - 1. )
    1040             !bb(j)=5.67e-8*skt(j)**4
    1041           end do
    1042 
    1043         do ibox=1,ncol
    1044           do j=1,npoints
    1045 
    1046             !add in surface emission
    1047 
    1048             fluxtop(j,ibox) = fluxtop(j,ibox)
    1049      &         + emsfc_lw * bb(j)
    1050      &         * trans_layers_above(j,ibox)
    1051            
    1052           end do
    1053         end do
    1054 
    1055 cIM
    1056 c        if (ncolprint.ne.0) then
    1057 
    1058 c          do j=1,npoints ,1000
    1059 c          write(6,'(a10)') 'j='
    1060 c          write(6,'(8I10)') j
    1061 c          write (6,'(a)') 'id:'
    1062 c          write (6,'(a)') 'surface'
    1063 
    1064 c          write (6,'(a)') 'emiss_layer:'
    1065 c          write (6,'(8f7.2)') (dem(1,ibox),ibox=1,ncolprint)
    1066    
    1067 c          write (6,'(a)') '100.*bb(j):'
    1068 c          write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint)
    1069    
    1070 c          write (6,'(a)') '100.*f:'
    1071 c          write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint)
    1072 c          end do
    1073 c        endif
    1074    
    1075         !now that you have the top of atmosphere radiance account
    1076         !for ISCCP procedures to determine cloud top temperature
    1077 
    1078         !account for partially transmitting cloud recompute flux
    1079         !ISCCP would see assuming a single layer cloud
    1080         !note choice here of 2.13, as it is primarily ice
    1081         !clouds which have partial emissivity and need the
    1082         !adjustment performed in this section
    1083         !
    1084         !If it turns out that the cloud brightness temperature
    1085         !is greater than 260K, then the liquid cloud conversion
    1086         !factor of 2.56 is used.
    1087         !
    1088         !Note that this is discussed on pages 85-87 of
    1089         !the ISCCP D level documentation (Rossow et al. 1996)
    1090            
    1091           do j=1,npoints 
    1092             !compute minimum brightness temperature and optical depth
    1093             btcmin(j) = 1. /  ( exp(1307.27/(attrop(j)-5.)) - 1. )
    1094           enddo
    1095         do ibox=1,ncol
    1096           do j=1,npoints 
    1097             transmax(j) = (fluxtop(j,ibox)-btcmin(j))
    1098      &                /(fluxtop_clrsky(j)-btcmin(j))
    1099             !note that the initial setting of tauir(j) is needed so that
    1100             !tauir(j) has a realistic value should the next if block be
    1101             !bypassed
    1102             tauir(j) = tau(j,ibox) * rec2p13
    1103             taumin(j) = -1. * log(max(min(transmax(j),0.9999999),0.001))
    1104 
    1105           enddo
    1106 
    1107           if (top_height .eq. 1) then
    1108             do j=1,npoints 
    1109               if (transmax(j) .gt. 0.001 .and.
    1110      &          transmax(j) .le. 0.9999999) then
    1111                 fluxtopinit(j) = fluxtop(j,ibox)
    1112                 tauir(j) = tau(j,ibox) *rec2p13
    1113               endif
    1114             enddo
    1115             do icycle=1,2
    1116               do j=1,npoints 
    1117                 if (tau(j,ibox) .gt. (tauchk            )) then
    1118                 if (transmax(j) .gt. 0.001 .and.
    1119      &            transmax(j) .le. 0.9999999) then
    1120                   emcld(j,ibox) = 1. - exp(-1. * tauir(j)  )
    1121                   fluxtop(j,ibox) = fluxtopinit(j) -   
    1122      &              ((1.-emcld(j,ibox))*fluxtop_clrsky(j))
    1123                   fluxtop(j,ibox)=max(1.E-06,
    1124      &              (fluxtop(j,ibox)/emcld(j,ibox)))
    1125                   tb(j,ibox)= 1307.27
    1126      &              / (log(1. + (1./fluxtop(j,ibox))))
    1127                   if (tb(j,ibox) .gt. 260.) then
    1128                     tauir(j) = tau(j,ibox) / 2.56
    1129                   end if                         
    1130                 end if
    1131                 end if
    1132               enddo
    1133             enddo
    1134                
    1135           endif
    1136        
    1137           do j=1,npoints
    1138             if (tau(j,ibox) .gt. (tauchk            )) then
    1139                 !cloudy box
    1140                 tb(j,ibox)= 1307.27/ (log(1. + (1./fluxtop(j,ibox))))
    1141                 if (top_height.eq.1.and.tauir(j).lt.taumin(j)) then
    1142                          tb(j,ibox) = attrop(j) - 5.
    1143                          tau(j,ibox) = 2.13*taumin(j)
    1144                 end if
    1145             else
    1146                 !clear sky brightness temperature
    1147                 tb(j,ibox) = 1307.27/(log(1.+(1./fluxtop_clrsky(j))))
    1148             end if
    1149           enddo ! j
    1150         enddo ! ibox
    1151 
    1152 cIM
    1153 c        if (ncolprint.ne.0) then
    1154 
    1155 c          do j=1,npoints,1000
    1156 c          write(6,'(a10)') 'j='
    1157 c          write(6,'(8I10)') j
    1158 
    1159 c          write (6,'(a)') 'attrop:'
    1160 c          write (6,'(8f7.2)') (attrop(j))
    1161    
    1162 c          write (6,'(a)') 'btcmin:'
    1163 c          write (6,'(8f7.2)') (btcmin(j))
    1164    
    1165 c          write (6,'(a)') 'fluxtop_clrsky*100:'
    1166 c          write (6,'(8f7.2)')
    1167 c     &      (100.*fluxtop_clrsky(j))
    1168 
    1169 c          write (6,'(a)') '100.*f_adj:'
    1170 c          write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint)
    1171    
    1172 c          write (6,'(a)') 'transmax:'
    1173 c          write (6,'(8f7.2)') (transmax(ibox),ibox=1,ncolprint)
    1174    
    1175 c          write (6,'(a)') 'tau:'
    1176 c          write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint)
    1177    
    1178 c          write (6,'(a)') 'emcld:'
    1179 c          write (6,'(8f7.2)') (emcld(j,ibox),ibox=1,ncolprint)
    1180    
    1181 c          write (6,'(a)') 'total_trans:'
    1182 c          write (6,'(8f7.2)')
    1183 c     &          (trans_layers_above(j,ibox),ibox=1,ncolprint)
    1184    
    1185 c          write (6,'(a)') 'total_emiss:'
    1186 c          write (6,'(8f7.2)')
    1187 c     &          (1.0-trans_layers_above(j,ibox),ibox=1,ncolprint)
    1188    
    1189 c          write (6,'(a)') 'total_trans:'
    1190 c          write (6,'(8f7.2)')
    1191 c     &          (trans_layers_above(j,ibox),ibox=1,ncolprint)
    1192    
    1193 c          write (6,'(a)') 'ppout:'
    1194 c          write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint)
    1195 c          enddo ! j
    1196 c        endif
    1197 
    1198       end if
    1199 
    1200 !     ---------------------------------------------------!
    1201 
    1202 !     
    1203 !     ---------------------------------------------------!
    1204 !     DETERMINE CLOUD TOP PRESSURE
    1205 !
    1206 !     again the 2 methods differ according to whether
    1207 !     or not you use the physical cloud top pressure (top_height = 2)
    1208 !     or the radiatively determined cloud top pressure (top_height = 1 or 3)
    1209 !
    1210 
    1211       !compute cloud top pressure
    1212       do 30 ibox=1,ncol
    1213         !segregate according to optical thickness
    1214         if (top_height .eq. 1 .or. top_height .eq. 3) then 
    1215           !find level whose temperature
    1216           !most closely matches brightness temperature
    1217           do j=1,npoints
    1218             nmatch(j)=0
    1219           enddo
    1220           do 29 ilev=1,nlev-1
    1221 !cdir nodep
    1222             do j=1,npoints
    1223               if ((at(j,ilev)   .ge. tb(j,ibox) .and.
    1224      &          at(j,ilev+1) .lt. tb(j,ibox)) .or.
    1225      &          (at(j,ilev) .le. tb(j,ibox) .and.
    1226      &          at(j,ilev+1) .gt. tb(j,ibox))) then
    1227    
    1228                 nmatch(j)=nmatch(j)+1
    1229                 if(abs(at(j,ilev)-tb(j,ibox)) .lt.
    1230      &            abs(at(j,ilev+1)-tb(j,ibox))) then
    1231                   match(j,nmatch(j))=ilev
    1232                 else
    1233                   match(j,nmatch(j))=ilev+1
    1234                 end if
    1235               end if                       
    1236             enddo
    1237 29        continue
    1238 
    1239           do j=1,npoints
    1240             if (nmatch(j) .ge. 1) then
    1241               ptop(j,ibox)=pfull(j,match(j,nmatch(j)))
    1242               levmatch(j,ibox)=match(j,nmatch(j))   
    1243             else
    1244               if (tb(j,ibox) .lt. atmin(j)) then
    1245                 ptop(j,ibox)=ptrop(j)
    1246                 levmatch(j,ibox)=itrop(j)
    1247               end if
    1248               if (tb(j,ibox) .gt. atmax(j)) then
    1249                 ptop(j,ibox)=pfull(j,nlev)
    1250                 levmatch(j,ibox)=nlev
    1251               end if                               
    1252             end if
    1253           enddo ! j
    1254 
    1255         else ! if (top_height .eq. 1 .or. top_height .eq. 3)
    1256  
    1257           do j=1,npoints     
    1258             ptop(j,ibox)=0.
    1259           enddo
    1260           do ilev=1,nlev
    1261             do j=1,npoints     
    1262               if ((ptop(j,ibox) .eq. 0. )
    1263 cIM  &           .and.(frac_out(j,ibox,ilev) .ne. 0)) then
    1264      &           .and.(frac_out(j,ibox,ilev) .ne. 0.0)) then
    1265                 ptop(j,ibox)=pfull(j,ilev)
    1266                 levmatch(j,ibox)=ilev
    1267               end if
    1268             end do
    1269           end do
    1270         end if                           
    1271          
    1272         do j=1,npoints
    1273           if (tau(j,ibox) .le. (tauchk            )) then
    1274             ptop(j,ibox)=0.
    1275             levmatch(j,ibox)=0     
    1276           endif
    1277         enddo
    1278 
    1279 30    continue
    1280              
    1281 !
    1282 !
    1283 !     ---------------------------------------------------!
    1284 
    1285 
    1286 !     
    1287 !     ---------------------------------------------------!
    1288 !     DETERMINE ISCCP CLOUD TYPE FREQUENCIES
    1289 !
    1290 !     Now that ptop and tau have been determined,
    1291 !     determine amount of each of the 49 ISCCP cloud
    1292 !     types
    1293 !
    1294 !     Also compute grid box mean cloud top pressure and
    1295 !     optical thickness.  The mean cloud top pressure and
    1296 !     optical thickness are averages over the cloudy
    1297 !     area only. The mean cloud top pressure is a linear
    1298 !     average of the cloud top pressures.  The mean cloud
    1299 !     optical thickness is computed by converting optical
    1300 !     thickness to an albedo, averaging in albedo units,
    1301 !     then converting the average albedo back to a mean
    1302 !     optical thickness. 
    1303 !
    1304 
    1305       !compute isccp frequencies
    1306 
    1307       !reset frequencies
    1308       do 38 ilev=1,7
    1309       do 38 ilev2=1,7
    1310         do j=1,npoints !
    1311              fq_isccp(j,ilev,ilev2)=0.
    1312         enddo
    1313 38    continue
    1314 
    1315       !reset variables need for averaging cloud properties
    1316       do j=1,npoints
    1317         totalcldarea(j) = 0.
    1318         meanalbedocld(j) = 0.
    1319         meanptop(j) = 0.
    1320         meantaucld(j) = 0.
    1321       enddo ! j
    1322 
    1323       boxarea = 1./real(ncol)
    1324      
    1325               !determine optical depth category
    1326 cIM       do 39 j=1,npoints
    1327 cIM       do ibox=1,ncol
    1328         do 39 ibox=1,ncol
    1329           do j=1,npoints
    1330 
    1331 cIM
    1332 c         CALL CPU_time(t1)
    1333 cIM
    1334 
    1335           if (tau(j,ibox) .gt. (tauchk            )
    1336      &      .and. ptop(j,ibox) .gt. 0.) then
    1337               box_cloudy(j,ibox)=.true.
    1338           endif
    1339 
    1340 cIM
    1341 c         CALL CPU_time(t2)
    1342 c         print*,'IF tau t2 - t1',t2 - t1
    1343 
    1344 c         CALL CPU_time(t1)
    1345 cIM
    1346 
    1347           if (box_cloudy(j,ibox)) then
    1348 
    1349               ! totalcldarea always diagnosed day or night
    1350               totalcldarea(j) = totalcldarea(j) + boxarea
    1351 
    1352               if (sunlit(j).eq.1) then
    1353 
    1354                 ! tau diagnostics only with sunlight
    1355 
    1356                 boxtau(j,ibox) = tau(j,ibox)
    1357 
    1358                 !convert optical thickness to albedo
    1359                   albedocld(j,ibox)
    1360      &            =real(invtau(min(nint(100.*tau(j,ibox)),45000)))
    1361            
    1362                 !contribute to averaging
    1363                 meanalbedocld(j) = meanalbedocld(j)
    1364      &            +albedocld(j,ibox)*boxarea
    1365 
    1366             endif
    1367 
    1368           endif
    1369          
    1370 cIM
    1371 c         CALL CPU_time(t2)
    1372 c         print*,'IF box_cloudy t2 - t1',t2 - t1
    1373          
    1374 c         CALL CPU_time(t1)
    1375 cIM BEG
    1376 cIM           !convert ptop to millibars
    1377               ptop(j,ibox)=ptop(j,ibox) / 100.
    1378            
    1379 cIM           !save for output cloud top pressure and optical thickness
    1380               boxptop(j,ibox) = ptop(j,ibox)
    1381 cIM END
    1382 
    1383 cIM BEG
    1384               !reset itau(j), ipres(j)
    1385               itau(j) = 0
    1386               ipres(j) = 0
    1387 
    1388               if (tau(j,ibox) .lt. isccp_taumin) then
    1389                   itau(j)=1
    1390               else if (tau(j,ibox) .ge. isccp_taumin
    1391      &                                   
    1392      &          .and. tau(j,ibox) .lt. 1.3) then
    1393                 itau(j)=2
    1394               else if (tau(j,ibox) .ge. 1.3
    1395      &          .and. tau(j,ibox) .lt. 3.6) then
    1396                 itau(j)=3
    1397               else if (tau(j,ibox) .ge. 3.6
    1398      &          .and. tau(j,ibox) .lt. 9.4) then
    1399                   itau(j)=4
    1400               else if (tau(j,ibox) .ge. 9.4
    1401      &          .and. tau(j,ibox) .lt. 23.) then
    1402                   itau(j)=5
    1403               else if (tau(j,ibox) .ge. 23.
    1404      &          .and. tau(j,ibox) .lt. 60.) then
    1405                   itau(j)=6
    1406               else if (tau(j,ibox) .ge. 60.) then
    1407                   itau(j)=7
    1408               end if
    1409 
    1410               !determine cloud top pressure category
    1411               if (    ptop(j,ibox) .gt. 0. 
    1412      &          .and.ptop(j,ibox) .lt. 180.) then
    1413                   ipres(j)=1
    1414               else if(ptop(j,ibox) .ge. 180.
    1415      &          .and.ptop(j,ibox) .lt. 310.) then
    1416                   ipres(j)=2
    1417               else if(ptop(j,ibox) .ge. 310.
    1418      &          .and.ptop(j,ibox) .lt. 440.) then
    1419                   ipres(j)=3
    1420               else if(ptop(j,ibox) .ge. 440.
    1421      &          .and.ptop(j,ibox) .lt. 560.) then
    1422                   ipres(j)=4
    1423               else if(ptop(j,ibox) .ge. 560.
    1424      &          .and.ptop(j,ibox) .lt. 680.) then
    1425                   ipres(j)=5
    1426               else if(ptop(j,ibox) .ge. 680.
    1427      &          .and.ptop(j,ibox) .lt. 800.) then
    1428                   ipres(j)=6
    1429               else if(ptop(j,ibox) .ge. 800.) then
    1430                   ipres(j)=7
    1431               end if
    1432 cIM END
    1433 
    1434           if (sunlit(j).eq.1 .or. top_height .eq. 3) then
    1435 
    1436 cIM         !convert ptop to millibars
    1437 cIM           ptop(j,ibox)=ptop(j,ibox) / 100.
    1438            
    1439 cIM         !save for output cloud top pressure and optical thickness
    1440 cIM             boxptop(j,ibox) = ptop(j,ibox)
    1441    
    1442             if (box_cloudy(j,ibox)) then
    1443            
    1444               meanptop(j) = meanptop(j) + ptop(j,ibox)*boxarea
    1445 
    1446 cIM             !reset itau(j), ipres(j)
    1447 cIM           itau(j) = 0
    1448 cIM           ipres(j) = 0
    1449 
    1450 c             if (tau(j,ibox) .lt. isccp_taumin) then
    1451 c                 itau(j)=1
    1452 c             else if (tau(j,ibox) .ge. isccp_taumin
    1453 c    &                                   
    1454 c    &          .and. tau(j,ibox) .lt. 1.3) then
    1455 c               itau(j)=2
    1456 c             else if (tau(j,ibox) .ge. 1.3
    1457 c    &          .and. tau(j,ibox) .lt. 3.6) then
    1458 c               itau(j)=3
    1459 c             else if (tau(j,ibox) .ge. 3.6
    1460 c    &          .and. tau(j,ibox) .lt. 9.4) then
    1461 c                 itau(j)=4
    1462 c             else if (tau(j,ibox) .ge. 9.4
    1463 c    &          .and. tau(j,ibox) .lt. 23.) then
    1464 c                 itau(j)=5
    1465 c             else if (tau(j,ibox) .ge. 23.
    1466 c    &          .and. tau(j,ibox) .lt. 60.) then
    1467 c                 itau(j)=6
    1468 c             else if (tau(j,ibox) .ge. 60.) then
    1469 c                 itau(j)=7
    1470 c             end if
    1471 
    1472 c             !determine cloud top pressure category
    1473 c             if (    ptop(j,ibox) .gt. 0. 
    1474 c    &          .and.ptop(j,ibox) .lt. 180.) then
    1475 c                 ipres(j)=1
    1476 c             else if(ptop(j,ibox) .ge. 180.
    1477 c    &          .and.ptop(j,ibox) .lt. 310.) then
    1478 c                 ipres(j)=2
    1479 c             else if(ptop(j,ibox) .ge. 310.
    1480 c    &          .and.ptop(j,ibox) .lt. 440.) then
    1481 c                 ipres(j)=3
    1482 c            else if(ptop(j,ibox) .ge. 440.
    1483 c    &          .and.ptop(j,ibox) .lt. 560.) then
    1484 c                 ipres(j)=4
    1485 c             else if(ptop(j,ibox) .ge. 560.
    1486 c    &          .and.ptop(j,ibox) .lt. 680.) then
    1487 c                 ipres(j)=5
    1488 c             else if(ptop(j,ibox) .ge. 680.
    1489 c    &          .and.ptop(j,ibox) .lt. 800.) then
    1490 c                 ipres(j)=6
    1491 c             else if(ptop(j,ibox) .ge. 800.) then
    1492 c                 ipres(j)=7
    1493 c             end if
    1494 
    1495               !update frequencies
    1496               if(ipres(j) .gt. 0.and.itau(j) .gt. 0) then
    1497               fq_isccp(j,itau(j),ipres(j))=
    1498      &          fq_isccp(j,itau(j),ipres(j))+ boxarea
    1499               end if
    1500 
    1501 cIM calcul stats regime dynamique BEG
    1502 !             iw(j) = int((w(j)-wmin)/pas_w) +1
    1503 !             pctj(itau(j),ipres(j),iw(j))=.FALSE.
    1504 !             !update frequencies W500
    1505 !             if (pct_ocean(j)) then
    1506 !             if (ipres(j) .gt. 0.and.itau(j) .gt. 0) then
    1507 !             if (iw(j) .gt. int(wmin).and.iw(j) .le. iwmx) then
    1508 c             print*,' ISCCP iw=',iw(j),j
    1509 !             fq_dynreg(itau(j),ipres(j),iw(j))=
    1510 !    &          fq_dynreg(itau(j),ipres(j),iw(j))+
    1511 !    &          boxarea
    1512 c    &          fq_isccp(j,itau(j),ipres(j))
    1513 !             pctj(itau(j),ipres(j),iw(j))=.TRUE.
    1514 c             nfq_dynreg(itau(j),ipres(j),iw(j))=
    1515 c    &          nfq_dynreg(itau(j),ipres(j),iw(j))+1.
    1516 !              end if
    1517 !             end if
    1518 !             end if
    1519 cIM calcul stats regime dynamique END
    1520             end if !IM boxcloudy
    1521 
    1522           end if !IM sunlit
    1523                        
    1524 cIM
    1525 c         CALL CPU_time(t2)
    1526 c         print*,'IF sunlit boxcloudy t2 - t1',t2 - t1
    1527 cIM
    1528         enddo !IM ibox/j
    1529 
    1530 
    1531 cIM ajout stats s/ W500 BEG
    1532 cIM ajout stats s/ W500 END
    1533 
    1534 c             if(j.EQ.6722) then
    1535 c             print*,' ISCCP',w(j),iw(j),ipres(j),itau(j)
    1536 c             endif
    1537 
    1538 !     if (pct_ocean(j)) then
    1539 !     if (ipres(j) .gt. 0.and.itau(j) .gt. 0) then
    1540 !     if (iw(j) .gt. int(wmin).and.iw(j) .le. iwmx) then
    1541 !     if(pctj(itau(j),ipres(j),iw(j))) THEN
    1542 !         nfq_dynreg(itau(j),ipres(j),iw(j))=
    1543 !    &    nfq_dynreg(itau(j),ipres(j),iw(j))+1.
    1544 c         if(itau(j).EQ.4.AND.ipres(j).EQ.2.AND.
    1545 c    &    iw(j).EQ.10) then
    1546 c         PRINT*,' isccp AVANT',
    1547 c    &    nfq_dynreg(itau(j),ipres(j),iw(j)),
    1548 c    &    fq_dynreg(itau(j),ipres(j),iw(j))
    1549 c         endif
    1550 !     endif
    1551 !     endif
    1552 !     endif
    1553 !     endif
    1554 39    continue !IM j/ibox
    1555      
    1556       !compute mean cloud properties
    1557       do j=1,npoints
    1558        if (totalcldarea(j) .gt. 0.) then
    1559          meanptop(j) = meanptop(j) / totalcldarea(j)
    1560          if (sunlit(j).eq.1) then
    1561            meanalbedocld(j)=meanalbedocld(j) / totalcldarea(j)
    1562            meantaucld(j)=tautab(min(255,max(1,nint(meanalbedocld(j)))))
    1563          end if
    1564        end if
    1565       enddo ! j
    1566 !
    1567 cIM ajout stats s/ W500 BEG
    1568 !     do nw = 1, iwmx
    1569 !     do l = 1, 7
    1570 !     do k = 1, 7
    1571 !       if (nfq_dynreg(k,l,nw).GT.0.) then
    1572 !       fq_dynreg(k,l,nw) = fq_dynreg(k,l,nw)/nfq_dynreg(k,l,nw)
    1573 c        if(k.EQ.4.AND.l.EQ.2.AND.nw.EQ.10) then
    1574 c        print*,' isccp APRES',nfq_dynreg(k,l,nw),
    1575 c    &   fq_dynreg(k,l,nw)
    1576 c        endif
    1577 !       else
    1578 !        if(fq_dynreg(k,l,nw).NE.0.) then
    1579 !        print*,'nfq_dynreg = 0 tau,pc,nw',k,l,nw,fq_dynreg(k,l,nw)
    1580 !        endif
    1581 c        fq_dynreg(k,l,nw) = -1.E+20
    1582 c        nfq_dynreg(k,l,nw) = 1.E+20
    1583 !       end if
    1584 !     enddo !k
    1585 !     enddo !l
    1586 !     enddo !nw
    1587 cIM ajout stats s/ W500 END
    1588 !     ---------------------------------------------------!
    1589 
    1590 !     ---------------------------------------------------!
    1591 !     OPTIONAL PRINTOUT OF DATA TO CHECK PROGRAM
    1592 !
    1593 !cIM
    1594 c      if (debugcol.ne.0) then
    1595 !     
    1596 c         do j=1,npoints,debugcol
    1597 
    1598 c            !produce character output
    1599 c            do ilev=1,nlev
    1600 c              do ibox=1,ncol
    1601 c                   acc(ilev,ibox)=0
    1602 c              enddo
    1603 c            enddo
    1604 
    1605 c            do ilev=1,nlev
    1606 c              do ibox=1,ncol
    1607 c                   acc(ilev,ibox)=frac_out(j,ibox,ilev)*2
    1608 c                   if (levmatch(j,ibox) .eq. ilev)
    1609 c     &                 acc(ilev,ibox)=acc(ilev,ibox)+1
    1610 c              enddo
    1611 c            enddo
    1612 
    1613              !print test
    1614 
    1615 c          write(ftn09,11) j
    1616 c11        format('ftn09.',i4.4)
    1617 c         open(9, FILE=ftn09, FORM='FORMATTED')
    1618 
    1619 c             write(9,'(a1)') ' '
    1620 c                    write(9,'(10i5)')
    1621 c     &                  (ilev,ilev=5,nlev,5)
    1622 c             write(9,'(a1)') ' '
    1623              
    1624 c             do ibox=1,ncol
    1625 c               write(9,'(40(a1),1x,40(a1))')
    1626 c     &           (cchar_realtops(acc(ilev,ibox)+1),ilev=1,nlev)
    1627 c     &           ,(cchar(acc(ilev,ibox)+1),ilev=1,nlev)
    1628 c             end do
    1629 c            close(9)
    1630 c
    1631 cIM
    1632 c             if (ncolprint.ne.0) then
    1633 c               write(6,'(a1)') ' '
    1634 c                    write(6,'(a2,1X,5(a7,1X),a50)')
    1635 c     &                  'ilev',
    1636 c     &                  'pfull','at',
    1637 c     &                  'cc*100','dem_s','dtau_s',
    1638 c     &                  'cchar'
    1639 
    1640 !               do 4012 ilev=1,nlev
    1641 !                    write(6,'(60i2)') (box(i,ilev),i=1,ncolprint)
    1642 !                   write(6,'(i2,1X,5(f7.2,1X),50(a1))')
    1643 !     &                  ilev,
    1644 !     &                  pfull(j,ilev)/100.,at(j,ilev),
    1645 !     &                  cc(j,ilev)*100.0,dem_s(j,ilev),dtau_s(j,ilev)
    1646 !     &                  ,(cchar(acc(ilev,ibox)+1),ibox=1,ncolprint)
    1647 !4012           continue
    1648 c               write (6,'(a)') 'skt(j):'
    1649 c               write (6,'(8f7.2)') skt(j)
    1650                                      
    1651 c               write (6,'(8I7)') (ibox,ibox=1,ncolprint)
    1652              
    1653 c               write (6,'(a)') 'tau:'
    1654 c               write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint)
    1655    
    1656 c               write (6,'(a)') 'tb:'
    1657 c               write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint)
    1658    
    1659 c               write (6,'(a)') 'ptop:'
    1660 c               write (6,'(8f7.2)') (ptop(j,ibox),ibox=1,ncolprint)
    1661 c             endif
    1662    
    1663 c        enddo
    1664        
    1665 c      end if
    1666 
    1667       return
    1668       end
     1185              match(j, nmatch(j)) = ilev + 1
     1186            END IF
     1187          END IF
     1188        END DO
     1189      END DO
     1190
     1191      DO j = 1, npoints
     1192        IF (nmatch(j)>=1) THEN
     1193          ptop(j, ibox) = pfull(j, match(j,nmatch(j)))
     1194          levmatch(j, ibox) = match(j, nmatch(j))
     1195        ELSE
     1196          IF (tb(j,ibox)<atmin(j)) THEN
     1197            ptop(j, ibox) = ptrop(j)
     1198            levmatch(j, ibox) = itrop(j)
     1199          END IF
     1200          IF (tb(j,ibox)>atmax(j)) THEN
     1201            ptop(j, ibox) = pfull(j, nlev)
     1202            levmatch(j, ibox) = nlev
     1203          END IF
     1204        END IF
     1205      END DO ! j
     1206
     1207    ELSE ! if (top_height .eq. 1 .or. top_height .eq. 3)
     1208
     1209      DO j = 1, npoints
     1210        ptop(j, ibox) = 0.
     1211      END DO
     1212      DO ilev = 1, nlev
     1213        DO j = 1, npoints
     1214          IF ((ptop(j,ibox)==0.) & ! IM  &
     1215                                   ! .and.(frac_out(j,ibox,ilev) .ne. 0))
     1216                                   ! then
     1217              .AND. (frac_out(j,ibox,ilev)/=0.0)) THEN
     1218            ptop(j, ibox) = pfull(j, ilev)
     1219            levmatch(j, ibox) = ilev
     1220          END IF
     1221        END DO
     1222      END DO
     1223    END IF
     1224
     1225    DO j = 1, npoints
     1226      IF (tau(j,ibox)<=(tauchk)) THEN
     1227        ptop(j, ibox) = 0.
     1228        levmatch(j, ibox) = 0
     1229      END IF
     1230    END DO
     1231
     1232  END DO
     1233
     1234
     1235
     1236  ! ---------------------------------------------------!
     1237
     1238
     1239
     1240  ! ---------------------------------------------------!
     1241  ! DETERMINE ISCCP CLOUD TYPE FREQUENCIES
     1242
     1243  ! Now that ptop and tau have been determined,
     1244  ! determine amount of each of the 49 ISCCP cloud
     1245  ! types
     1246
     1247  ! Also compute grid box mean cloud top pressure and
     1248  ! optical thickness.  The mean cloud top pressure and
     1249  ! optical thickness are averages over the cloudy
     1250  ! area only. The mean cloud top pressure is a linear
     1251  ! average of the cloud top pressures.  The mean cloud
     1252  ! optical thickness is computed by converting optical
     1253  ! thickness to an albedo, averaging in albedo units,
     1254  ! then converting the average albedo back to a mean
     1255  ! optical thickness.
     1256
     1257
     1258    !compute isccp frequencies
     1259
     1260    !reset frequencies
     1261  DO ilev = 1, 7
     1262    DO ilev2 = 1, 7
     1263      DO j = 1, npoints !
     1264        fq_isccp(j, ilev, ilev2) = 0.
     1265      END DO
     1266    END DO
     1267  END DO
     1268
     1269    !reset variables need for averaging cloud properties
     1270  DO j = 1, npoints
     1271    totalcldarea(j) = 0.
     1272    meanalbedocld(j) = 0.
     1273    meanptop(j) = 0.
     1274    meantaucld(j) = 0.
     1275  END DO ! j
     1276
     1277  boxarea = 1./real(ncol)
     1278
     1279    !determine optical depth category
     1280  ! IM       do 39 j=1,npoints
     1281  ! IM       do ibox=1,ncol
     1282  DO ibox = 1, ncol
     1283    DO j = 1, npoints
     1284
     1285      ! IM
     1286      ! CALL CPU_time(t1)
     1287      ! IM
     1288
     1289      IF (tau(j,ibox)>(tauchk) .AND. ptop(j,ibox)>0.) THEN
     1290        box_cloudy(j, ibox) = .TRUE.
     1291      END IF
     1292
     1293      ! IM
     1294      ! CALL CPU_time(t2)
     1295      ! print*,'IF tau t2 - t1',t2 - t1
     1296
     1297      ! CALL CPU_time(t1)
     1298      ! IM
     1299
     1300      IF (box_cloudy(j,ibox)) THEN
     1301
     1302          ! totalcldarea always diagnosed day or night
     1303        totalcldarea(j) = totalcldarea(j) + boxarea
     1304
     1305        IF (sunlit(j)==1) THEN
     1306
     1307            ! tau diagnostics only with sunlight
     1308
     1309          boxtau(j, ibox) = tau(j, ibox)
     1310
     1311            !convert optical thickness to albedo
     1312          albedocld(j, ibox) = real(invtau(min(nint(100.*tau(j,ibox)), &
     1313            45000)))
     1314
     1315            !contribute to averaging
     1316          meanalbedocld(j) = meanalbedocld(j) + albedocld(j, ibox)*boxarea
     1317
     1318        END IF
     1319
     1320      END IF
     1321
     1322      ! IM
     1323      ! CALL CPU_time(t2)
     1324      ! print*,'IF box_cloudy t2 - t1',t2 - t1
     1325
     1326      ! CALL CPU_time(t1)
     1327      ! IM BEG
     1328      ! IM           !convert ptop to millibars
     1329      ptop(j, ibox) = ptop(j, ibox)/100.
     1330
     1331      ! IM           !save for output cloud top pressure and optical
     1332      ! thickness
     1333      boxptop(j, ibox) = ptop(j, ibox)
     1334      ! IM END
     1335
     1336      ! IM BEG
     1337        !reset itau(j), ipres(j)
     1338      itau(j) = 0
     1339      ipres(j) = 0
     1340
     1341      IF (tau(j,ibox)<isccp_taumin) THEN
     1342        itau(j) = 1
     1343      ELSE IF (tau(j,ibox)>=isccp_taumin .AND. tau(j,ibox)<1.3) THEN
     1344        itau(j) = 2
     1345      ELSE IF (tau(j,ibox)>=1.3 .AND. tau(j,ibox)<3.6) THEN
     1346        itau(j) = 3
     1347      ELSE IF (tau(j,ibox)>=3.6 .AND. tau(j,ibox)<9.4) THEN
     1348        itau(j) = 4
     1349      ELSE IF (tau(j,ibox)>=9.4 .AND. tau(j,ibox)<23.) THEN
     1350        itau(j) = 5
     1351      ELSE IF (tau(j,ibox)>=23. .AND. tau(j,ibox)<60.) THEN
     1352        itau(j) = 6
     1353      ELSE IF (tau(j,ibox)>=60.) THEN
     1354        itau(j) = 7
     1355      END IF
     1356
     1357        !determine cloud top pressure category
     1358      IF (ptop(j,ibox)>0. .AND. ptop(j,ibox)<180.) THEN
     1359        ipres(j) = 1
     1360      ELSE IF (ptop(j,ibox)>=180. .AND. ptop(j,ibox)<310.) THEN
     1361        ipres(j) = 2
     1362      ELSE IF (ptop(j,ibox)>=310. .AND. ptop(j,ibox)<440.) THEN
     1363        ipres(j) = 3
     1364      ELSE IF (ptop(j,ibox)>=440. .AND. ptop(j,ibox)<560.) THEN
     1365        ipres(j) = 4
     1366      ELSE IF (ptop(j,ibox)>=560. .AND. ptop(j,ibox)<680.) THEN
     1367        ipres(j) = 5
     1368      ELSE IF (ptop(j,ibox)>=680. .AND. ptop(j,ibox)<800.) THEN
     1369        ipres(j) = 6
     1370      ELSE IF (ptop(j,ibox)>=800.) THEN
     1371        ipres(j) = 7
     1372      END IF
     1373      ! IM END
     1374
     1375      IF (sunlit(j)==1 .OR. top_height==3) THEN
     1376
     1377        ! IM         !convert ptop to millibars
     1378        ! IM           ptop(j,ibox)=ptop(j,ibox) / 100.
     1379
     1380        ! IM         !save for output cloud top pressure and optical
     1381        ! thickness
     1382        ! IM             boxptop(j,ibox) = ptop(j,ibox)
     1383
     1384        IF (box_cloudy(j,ibox)) THEN
     1385
     1386          meanptop(j) = meanptop(j) + ptop(j, ibox)*boxarea
     1387
     1388          ! IM             !reset itau(j), ipres(j)
     1389          ! IM           itau(j) = 0
     1390          ! IM           ipres(j) = 0
     1391
     1392          ! if (tau(j,ibox) .lt. isccp_taumin) then
     1393          ! itau(j)=1
     1394          ! else if (tau(j,ibox) .ge. isccp_taumin
     1395          ! &
     1396          ! &          .and. tau(j,ibox) .lt. 1.3) then
     1397          ! itau(j)=2
     1398          ! else if (tau(j,ibox) .ge. 1.3
     1399          ! &          .and. tau(j,ibox) .lt. 3.6) then
     1400          ! itau(j)=3
     1401          ! else if (tau(j,ibox) .ge. 3.6
     1402          ! &          .and. tau(j,ibox) .lt. 9.4) then
     1403          ! itau(j)=4
     1404          ! else if (tau(j,ibox) .ge. 9.4
     1405          ! &          .and. tau(j,ibox) .lt. 23.) then
     1406          ! itau(j)=5
     1407          ! else if (tau(j,ibox) .ge. 23.
     1408          ! &          .and. tau(j,ibox) .lt. 60.) then
     1409          ! itau(j)=6
     1410          ! else if (tau(j,ibox) .ge. 60.) then
     1411          ! itau(j)=7
     1412          ! end if
     1413
     1414          ! !determine cloud top pressure category
     1415          ! if (    ptop(j,ibox) .gt. 0.
     1416          ! &          .and.ptop(j,ibox) .lt. 180.) then
     1417          ! ipres(j)=1
     1418          ! else if(ptop(j,ibox) .ge. 180.
     1419          ! &          .and.ptop(j,ibox) .lt. 310.) then
     1420          ! ipres(j)=2
     1421          ! else if(ptop(j,ibox) .ge. 310.
     1422          ! &          .and.ptop(j,ibox) .lt. 440.) then
     1423          ! ipres(j)=3
     1424          ! else if(ptop(j,ibox) .ge. 440.
     1425          ! &          .and.ptop(j,ibox) .lt. 560.) then
     1426          ! ipres(j)=4
     1427          ! else if(ptop(j,ibox) .ge. 560.
     1428          ! &          .and.ptop(j,ibox) .lt. 680.) then
     1429          ! ipres(j)=5
     1430          ! else if(ptop(j,ibox) .ge. 680.
     1431          ! &          .and.ptop(j,ibox) .lt. 800.) then
     1432          ! ipres(j)=6
     1433          ! else if(ptop(j,ibox) .ge. 800.) then
     1434          ! ipres(j)=7
     1435          ! end if
     1436
     1437            !update frequencies
     1438          IF (ipres(j)>0 .AND. itau(j)>0) THEN
     1439            fq_isccp(j, itau(j), ipres(j)) = fq_isccp(j, itau(j), ipres(j)) + &
     1440              boxarea
     1441          END IF
     1442
     1443          ! IM calcul stats regime dynamique BEG
     1444          ! iw(j) = int((w(j)-wmin)/pas_w) +1
     1445          ! pctj(itau(j),ipres(j),iw(j))=.FALSE.
     1446          ! !update frequencies W500
     1447          ! if (pct_ocean(j)) then
     1448          ! if (ipres(j) .gt. 0.and.itau(j) .gt. 0) then
     1449          ! if (iw(j) .gt. int(wmin).and.iw(j) .le. iwmx) then
     1450          ! print*,' ISCCP iw=',iw(j),j
     1451          ! fq_dynreg(itau(j),ipres(j),iw(j))=
     1452          ! &          fq_dynreg(itau(j),ipres(j),iw(j))+
     1453          ! &          boxarea
     1454          ! &          fq_isccp(j,itau(j),ipres(j))
     1455          ! pctj(itau(j),ipres(j),iw(j))=.TRUE.
     1456          ! nfq_dynreg(itau(j),ipres(j),iw(j))=
     1457          ! &          nfq_dynreg(itau(j),ipres(j),iw(j))+1.
     1458          ! end if
     1459          ! end if
     1460          ! end if
     1461          ! IM calcul stats regime dynamique END
     1462        END IF !IM boxcloudy
     1463
     1464      END IF !IM sunlit
     1465
     1466      ! IM
     1467      ! CALL CPU_time(t2)
     1468      ! print*,'IF sunlit boxcloudy t2 - t1',t2 - t1
     1469      ! IM
     1470    END DO !IM ibox/j
     1471
     1472
     1473    ! IM ajout stats s/ W500 BEG
     1474    ! IM ajout stats s/ W500 END
     1475
     1476    ! if(j.EQ.6722) then
     1477    ! print*,' ISCCP',w(j),iw(j),ipres(j),itau(j)
     1478    ! endif
     1479
     1480    ! if (pct_ocean(j)) then
     1481    ! if (ipres(j) .gt. 0.and.itau(j) .gt. 0) then
     1482    ! if (iw(j) .gt. int(wmin).and.iw(j) .le. iwmx) then
     1483    ! if(pctj(itau(j),ipres(j),iw(j))) THEN
     1484    ! nfq_dynreg(itau(j),ipres(j),iw(j))=
     1485    ! &    nfq_dynreg(itau(j),ipres(j),iw(j))+1.
     1486    ! if(itau(j).EQ.4.AND.ipres(j).EQ.2.AND.
     1487    ! &    iw(j).EQ.10) then
     1488    ! PRINT*,' isccp AVANT',
     1489    ! &    nfq_dynreg(itau(j),ipres(j),iw(j)),
     1490    ! &    fq_dynreg(itau(j),ipres(j),iw(j))
     1491    ! endif
     1492    ! endif
     1493    ! endif
     1494    ! endif
     1495    ! endif
     1496
     1497  END DO
     1498    !compute mean cloud properties
     1499  ! IM j/ibox
     1500  DO j = 1, npoints
     1501    IF (totalcldarea(j)>0.) THEN
     1502      meanptop(j) = meanptop(j)/totalcldarea(j)
     1503      IF (sunlit(j)==1) THEN
     1504        meanalbedocld(j) = meanalbedocld(j)/totalcldarea(j)
     1505        meantaucld(j) = tautab(min(255,max(1,nint(meanalbedocld(j)))))
     1506      END IF
     1507    END IF
     1508  END DO ! j
     1509
     1510  ! IM ajout stats s/ W500 BEG
     1511  ! do nw = 1, iwmx
     1512  ! do l = 1, 7
     1513  ! do k = 1, 7
     1514  ! if (nfq_dynreg(k,l,nw).GT.0.) then
     1515  ! fq_dynreg(k,l,nw) = fq_dynreg(k,l,nw)/nfq_dynreg(k,l,nw)
     1516  ! if(k.EQ.4.AND.l.EQ.2.AND.nw.EQ.10) then
     1517  ! print*,' isccp APRES',nfq_dynreg(k,l,nw),
     1518  ! &   fq_dynreg(k,l,nw)
     1519  ! endif
     1520  ! else
     1521  ! if(fq_dynreg(k,l,nw).NE.0.) then
     1522  ! print*,'nfq_dynreg = 0 tau,pc,nw',k,l,nw,fq_dynreg(k,l,nw)
     1523  ! endif
     1524  ! fq_dynreg(k,l,nw) = -1.E+20
     1525  ! nfq_dynreg(k,l,nw) = 1.E+20
     1526  ! end if
     1527  ! enddo !k
     1528  ! enddo !l
     1529  ! enddo !nw
     1530  ! IM ajout stats s/ W500 END
     1531  ! ---------------------------------------------------!
     1532
     1533  ! ---------------------------------------------------!
     1534  ! OPTIONAL PRINTOUT OF DATA TO CHECK PROGRAM
     1535
     1536  ! cIM
     1537  ! if (debugcol.ne.0) then
     1538
     1539  ! do j=1,npoints,debugcol
     1540
     1541  ! !produce character output
     1542  ! do ilev=1,nlev
     1543  ! do ibox=1,ncol
     1544  ! acc(ilev,ibox)=0
     1545  ! enddo
     1546  ! enddo
     1547
     1548  ! do ilev=1,nlev
     1549  ! do ibox=1,ncol
     1550  ! acc(ilev,ibox)=frac_out(j,ibox,ilev)*2
     1551  ! if (levmatch(j,ibox) .eq. ilev)
     1552  ! &                 acc(ilev,ibox)=acc(ilev,ibox)+1
     1553  ! enddo
     1554  ! enddo
     1555
     1556    !print test
     1557
     1558  ! write(ftn09,11) j
     1559  ! 11        format('ftn09.',i4.4)
     1560  ! open(9, FILE=ftn09, FORM='FORMATTED')
     1561
     1562  ! write(9,'(a1)') ' '
     1563  ! write(9,'(10i5)')
     1564  ! &                  (ilev,ilev=5,nlev,5)
     1565  ! write(9,'(a1)') ' '
     1566
     1567  ! do ibox=1,ncol
     1568  ! write(9,'(40(a1),1x,40(a1))')
     1569  ! &           (cchar_realtops(acc(ilev,ibox)+1),ilev=1,nlev)
     1570  ! &           ,(cchar(acc(ilev,ibox)+1),ilev=1,nlev)
     1571  ! end do
     1572  ! close(9)
     1573
     1574  ! IM
     1575  ! if (ncolprint.ne.0) then
     1576  ! write(6,'(a1)') ' '
     1577  ! write(6,'(a2,1X,5(a7,1X),a50)')
     1578  ! &                  'ilev',
     1579  ! &                  'pfull','at',
     1580  ! &                  'cc*100','dem_s','dtau_s',
     1581  ! &                  'cchar'
     1582
     1583  ! do 4012 ilev=1,nlev
     1584  ! write(6,'(60i2)') (box(i,ilev),i=1,ncolprint)
     1585  ! write(6,'(i2,1X,5(f7.2,1X),50(a1))')
     1586  ! &                  ilev,
     1587  ! &                  pfull(j,ilev)/100.,at(j,ilev),
     1588  ! &                  cc(j,ilev)*100.0,dem_s(j,ilev),dtau_s(j,ilev)
     1589  ! &                  ,(cchar(acc(ilev,ibox)+1),ibox=1,ncolprint)
     1590  ! 4012           continue
     1591  ! write (6,'(a)') 'skt(j):'
     1592  ! write (6,'(8f7.2)') skt(j)
     1593
     1594  ! write (6,'(8I7)') (ibox,ibox=1,ncolprint)
     1595
     1596  ! write (6,'(a)') 'tau:'
     1597  ! write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint)
     1598
     1599  ! write (6,'(a)') 'tb:'
     1600  ! write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint)
     1601
     1602  ! write (6,'(a)') 'ptop:'
     1603  ! write (6,'(8f7.2)') (ptop(j,ibox),ibox=1,ncolprint)
     1604  ! endif
     1605
     1606  ! enddo
     1607
     1608  ! end if
     1609
     1610  RETURN
     1611END SUBROUTINE isccp_cloud_types
Note: See TracChangeset for help on using the changeset viewer.