concld(i,k) = min(rh(i,k),min(1._r8,max(0._r8,zdu(i,k)*5.e4)))         endif      end do   end do#endif!! Evaluate effective column-integrated convective cloud cover using! random overlap assumption (for diagnostic purposes only)!   do i=1,ncol      clrsky(i) = 1.0   end do   do k=pver,1,-1      do i=1,ncol         clrsky(i) = clrsky(i)*(1. - concld(i,k))      end do   end do   do i=1,ncol      clc(i) = 1. - clrsky(i)   end do!!          ****** Compute layer cloudiness ******!! There is effecively no top for high cloud formation (can for all the way! up to 1mb)! The bottom of middle level cloud (or the top of low level cloud) is! arbitrarily define to be 750 mb (premib)!   premib = 750.e2   pretop = 1.0e2                 ! top of cloud layer is at 1 mb!! Find most stable level below 750 mb for evaluating stratus regimes!   do i=1,ncol      dtdpmn(i) = 0.0      kdthdp(i) = 0      dthtdp(i,1) = 0.0   end do   do k=2,pver-2      do i=1,ncol         if (pmid(i,k) >= premib) then            dthdp = 100.0*(theta(i,k) - theta(i,k-1))*rpdeli(i,k-1)         else            dthdp = 0.0         end if         if (dthdp < dtdpmn(i)) then            dtdpmn(i) = dthdp            kdthdp(i) = k     ! index of interface of max inversion         end if         dthtdp(i,k) = dthdp      end do   end do   do k=pver-1,pver      do i=1,ncol         if (0.0 < dtdpmn(i)) then            dtdpmn(i) = 0.0         end if         dthtdp(i,k) = 0.0      end do   end do!! For some reason, allowing clouds in the bottom model layer causes bad! error growth characteristics!#ifdef PERGRO   numkcld = pver - 1#else   numkcld = pver#endif!! bvf => brunt-vaisalla frequency (approx. 1-sided diff.)! this stability measure is used to set a local relative humidity! threshold when evaluating the fractional area of layered cloud!   do 10 k=2,numkcld      kp1 = min(k + 1,pver)      do i=1,ncol         if (dthtdp(i,k) > dtdpmn(i)) then            dthtdp(i,k) = 0.         end if         cldbnd(i) = pmid(i,k).ge.pretop         lol(i) = pmid(i,k).ge.premib         rho = pmid(i,k)/(rair*temp(i,k))         bvf = -rho*gravit*gravit*((theta(i,k)-theta(i,k-1))* &                  rpdeli(i,k-1))/theta(i,k)         if (cldbnd(i)) then            rhlim = 0.999 - (1.0-rhminh)*(1.0-min(1.0_r8,max(0.0_r8,bvf*rbvflim)))            rhden = 1.0 - rhlim         else            rhlim = 0.999            rhden = 0.001         end if         rhdif = (rh(i,k) - rhlim)/rhden#ifdef PERGRO         cld9(i) = 0.1#else         cld9(i) = min(0.999_r8,(max(rhdif,0.0_r8))**2)#endif!! Ignore brunt-vaisalla stability estimate of local relative humidity! threshold when evaluating low cloud where local vertical motion is! less than some prescribed value (see low cloud section below)! Relative humidity threshold is fixed at rhminl for this case, except! over snow-free land, where it is reduced by 10%.  This distinction is! made to account for enhanced cloud drop nucleation ({\it i.e.,} at! lower relative humidities) that can occur over CCN rich land areas.!         if (lol(i)) then            if (land(i) .and. (snowh(i) <= 0.000001)) then               rhlim = rhminl - 0.10            else               rhlim = rhminl            endif            rhdif2 = (rh(i,k) - rhlim)/(1.0-rhlim)            cld8(i) = min(0.999_r8,(max(rhdif2,0.0_r8))**2)         else            cld8(i) = cld9(i)         end if!! save rhlim to rhu00, it handles well by itself for low/high cloud!         rhu00(i,k)=rhlim      end do!! Final evaluation of layered cloud fraction!      do i=1,ncol!! Low cloud: non-zero only if vertical velocity requirements are satisfied! Current vertical velocity threshold is omega < +50 mb/day with a 50 mb/day! linear ramp (other quantities in the class of "disposable" parameters)!         if (lol(i)) then            if (omga(i,k) < 0.05787) then               cld = cld8(i)* min(1.0_r8,max(0.0_r8,(0.05787-omga(i,k))/0.05787))            else               cld = 0.0!! give a fake value of rhlim, 2.0, which would never be used!                           rhu00(i,k)=2.0            end if            cloud(i,k) = cld#ifdef OLDLOWCLD!! Compute cloud associated with low level inversions.!            strat = max(0.,min(0.95_r8,-6.67*dthtdp(i,k) - 0.667))            rhb   = 1.0 - (0.9 - rh(i,k+1))/0.3            if (rh(i,k+1) < 0.6) then               strat = 0.0            end if            if (rh(i,k+1) >= 0.6 .and. rh(i,k+1) <= 0.9) then               strat = strat*rhb            end if!! Linear transition from stratus to trade cu as inversion rises.! Transition starts at 900mb and completes at 750mb (premib)!            pdepth = max(pmid(i,k) - 750.e2,0.0)            stratfac = min(pdepth,150.0e2_r8)/150.e2            if (dthtdp(i,k) <= -0.125 ) then               cloud(i,k) = strat*stratfac            else               cloud(i,k) = cld            end if#endif         else                  ! Middle and high level cloud            if ( cldbnd(i) ) then               cloud(i,k) = cld9(i)            else               cloud(i,k) = 0.0!! set a fake value for rhlim, 2.0!               rhu00(i,k)=2.0            end if         end if      end do10    continue                  ! k=2,pver-1#ifdef PERGRO      rhu00(:ncol,pver)=0.0#endif!! Add in the marine strat!      do i=1,ncol!!jrbee bugfix?!         if (kdthdp(i) /= 0) then            k = kdthdp(i)            kp1 = min(k+1,pver)            strat = min(1._r8,max(0._r8,(theta(i,k700)-thetas(i))*.057-.5573))!! assign the stratus to the layer just below max inversion! the relative humidity changes so rapidly across the inversion! that it is not safe to just look immediately below the inversion! so limit the stratus cloud by rh in both layers below the inversion!            if (ocean(i) .and. dthtdp(i,k) <= -0.125 ) then               cldst(i,k) = min(strat,max(rh(i,k),rh(i,kp1)))               cloud(i,k) = max(cloud(i,k),cldst(i,k))            else               cldst(i,k) = 0.            end if         end if      end do!! Merge convective and layered cloud fraction for total cloud!      do k=1,pver         do i=1,ncol!!          cloud(i,k) = max(0.0,min(0.999_r8,!     $                 (1.0 - concld(i,k))*cloud(i,k) + concld(i,k)))! change to a max overlap assumption between convective and strat clouds!            cloud(i,k) = max(0.0_r8,min(0.999_r8,max(concld(i,k),cloud(i,k))))#ifndef PERGRO            if (rh(i,k) > 0.99) then               cloud(i,k) = max(0.01_r8,cloud(i,k))            endif#endif         end do      end do!      returnend subroutine cldfrc