LCOV - code coverage report
Current view: | top level - physics/rrtmgp - mcica_subcol_gen.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 83 | 83 | 100.0 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 2 | 2 | 100.0 %

              Line data    Source code
           1             : module mcica_subcol_gen
           2             :
           3             : !----------------------------------------------------------------------------------------
           4             : !
           5             : ! Purpose: Create McICA stochastic arrays for cloud optical properties.
           6             : ! Input cloud optical properties directly: cloud optical depth, single
           7             : ! scattering albedo and asymmetry parameter.  Output will be stochastic
           8             : ! arrays of these variables.  (longwave scattering is not yet available)
           9             : !
          10             : ! Original code: From RRTMG, with the following copyright notice,
          11             : ! based on Raisanen et al., QJRMS, 2004:
          12             : !  --------------------------------------------------------------------------
          13             : ! |                                                                          |
          14             : ! |  Copyright 2006-2007, Atmospheric & Environmental Research, Inc. (AER).  |
          15             : ! |  This software may be used, copied, or redistributed as long as it is    |
          16             : ! |  not sold and this copyright notice is reproduced on each copy made.     |
          17             : ! |  This model is provided as is without any express or implied warranties. |
          18             : ! |                       (http://www.rtweb.aer.com/)                        |
          19             : ! |                                                                          |
          20             : !  --------------------------------------------------------------------------
          21             : ! This code is a refactored version of code originally in the files
          22             : ! mcica_subcol_gen_lw.F90 and mcica_subcol_gen_sw.F90
          23             : !
          24             : ! Uses the KISS random number generator.
          25             : !
          26             : ! Overlap assumption: maximum-random.
          27             : !
          28             : !----------------------------------------------------------------------------------------
          29             :
          30             : use shr_kind_mod,         only: r8 => shr_kind_r8
          31             : use ppgrid,               only: pcols, pver
          32             : use shr_RandNum_mod,      only: ShrKissRandGen
          33             : use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp
          34             :
          35             : implicit none
          36             : private
          37             : save
          38             :
          39             : public :: mcica_subcol_lw, mcica_subcol_sw
          40             :
          41             : !========================================================================================
          42             : contains
          43             : !========================================================================================
          44             :
          45      746136 : subroutine mcica_subcol_lw( &
          46             :    kdist, nbnd, ngpt, ncol, nver,           &
          47      746136 :    changeseed, pmid, cldfrac, tauc, taucmcl )
          48             :
          49             :    ! Arrays use CAM vertical index convention: index increases from top to bottom.
          50             :    ! This index ordering is assumed in the maximum-random overlap algorithm which starts
          51             :    ! at the top of a column and marches down, with each layer depending on the state
          52             :    ! of the layer above it.
          53             :    !
          54             :    ! For GCM mode, changeseed must be offset between LW and SW by at least the
          55             :    ! number of subcolumns
          56             :
          57             :    ! arguments
          58             :    class(ty_gas_optics_rrtmgp), intent(in) :: kdist  ! spectral information
          59             :    integer,  intent(in)  :: nbnd                     ! number of spectral bands
          60             :    integer,  intent(in)  :: ngpt                     ! number of subcolumns (g-point intervals)
          61             :    integer,  intent(in)  :: ncol                     ! number of columns
          62             :    integer,  intent(in)  :: nver                     ! number of layers
          63             :    integer,  intent(in)  :: changeseed               ! if the subcolumn generator is called multiple times,
          64             :                                                      ! permute the seed between each call.
          65             :    real(r8), intent(in)  :: pmid(pcols,pver)         ! layer pressures (Pa)
          66             :    real(r8), intent(in)  :: cldfrac(ncol,nver)       ! layer cloud fraction
          67             :    real(r8), intent(in)  :: tauc(nbnd,ncol,nver)     ! cloud optical depth
          68             :    real(r8), intent(out) :: taucmcl(ngpt,ncol,nver)  ! subcolumn cloud optical depth [mcica]
          69             :
          70             :    ! Local variables
          71             :
          72             :    integer :: i, isubcol, k, n
          73             :
          74             :    real(r8), parameter :: cldmin = 1.0e-80_r8  ! min cloud fraction
          75     1492272 :    real(r8) :: cldf(ncol,pver)      ! cloud fraction clipped to cldmin
          76             :
          77      746136 :    type(ShrKissRandGen) :: kiss_gen  ! KISS RNG object
          78     1492272 :    integer  :: kiss_seed(ncol,4)
          79     1492272 :    real(r8) :: rand_num_1d(ncol,1)   ! random number (kissvec)
          80     1492272 :    real(r8) :: rand_num(ncol,nver)   ! random number (kissvec)
          81             :
          82     1492272 :    real(r8) :: cdf(ngpt,ncol,nver)   ! random numbers
          83     1492272 :    logical  :: iscloudy(ngpt,ncol,nver)   ! flag that says whether a gridbox is cloudy
          84             :    !------------------------------------------------------------------------------------------
          85             :
          86             :    ! clip cloud fraction
          87  1159408584 :    cldf(:,:) = cldfrac(:ncol,:)
          88  1159408584 :    where (cldf(:,:) < cldmin)
          89             :       cldf(:,:) = 0._r8
          90             :    end where
          91             :
          92             :    ! Create a seed that depends on the state of the columns.
          93             :    ! Use pmid from bottom four layers.
          94    12458736 :    do i = 1, ncol
          95    11712600 :       kiss_seed(i,1) = (pmid(i,pver)   - int(pmid(i,pver)))    * 1000000000
          96    11712600 :       kiss_seed(i,2) = (pmid(i,pver-1) - int(pmid(i,pver-1)))  * 1000000000
          97    11712600 :       kiss_seed(i,3) = (pmid(i,pver-2) - int(pmid(i,pver-2)))  * 1000000000
          98    12458736 :       kiss_seed(i,4) = (pmid(i,pver-3) - int(pmid(i,pver-3)))  * 1000000000
          99             :    end do
         100             :
         101             :    ! create the RNG object
         102      746136 :    kiss_gen = ShrKissRandGen(kiss_seed)
         103             :
         104             :    ! Advance randum number generator by changeseed values
         105    96251544 :    do i = 1, changeSeed
         106    96251544 :       call kiss_gen%random(rand_num_1d)
         107             :    end do
         108             :
         109             :    ! Generate random numbers in each subcolumn at every level
         110    96251544 :    do isubcol = 1,ngpt
         111    95505408 :       call kiss_gen%random(rand_num)
         112 >14840*10^7 :       cdf(isubcol,:,:) = rand_num(:,:)
         113             :    enddo
         114             :
         115             :    ! Maximum-Random overlap
         116             :    ! i) pick a random number for top layer.
         117             :    ! ii) walk down the column:
         118             :    !    - if the layer above is cloudy, use the same random number as in the layer above
         119             :    !    - if the layer above is clear, use a new random number
         120             :
         121    69390648 :    do k = 2, nver
         122  1146949848 :       do i = 1, ncol
         123 >13907*10^7 :          do isubcol = 1, ngpt
         124 >13900*10^7 :             if (cdf(isubcol,i,k-1) > 1._r8 - cldf(i,k-1) ) then
         125 11091485592 :                cdf(isubcol,i,k) = cdf(isubcol,i,k-1)
         126             :             else
         127 >12683*10^7 :                cdf(isubcol,i,k) = cdf(isubcol,i,k) * (1._r8 - cldf(i,k-1))
         128             :             end if
         129             :          end do
         130             :       end do
         131             :    end do
         132             :
         133    70136784 :    do k = 1, nver
         134 >14058*10^7 :       iscloudy(:,:,k) = (cdf(:,:,k) >= 1._r8 - spread(cldf(:,k), dim=1, nCopies=ngpt) )
         135             :    end do
         136             :
         137             :    ! -- generate subcolumns for homogeneous clouds -----
         138             :    ! where there is a cloud, set the subcolumn cloud properties;
         139             :    ! incoming tauc should be in-cloud quantites and not grid-averaged quantities
         140    70136784 :    do k = 1,nver
         141  1159408584 :       do i = 1,ncol
         142 >14058*10^7 :          do isubcol = 1,ngpt
         143 >14051*10^7 :             if (iscloudy(isubcol,i,k) .and. (cldf(i,k) > 0._r8) ) then
         144 11261941296 :                n = kdist%convert_gpt2band(isubcol)
         145 11261941296 :                taucmcl(isubcol,i,k) = tauc(n,i,k)
         146             :             else
         147 >12816*10^7 :                taucmcl(isubcol,i,k) = 0._r8
         148             :             end if
         149             :          end do
         150             :       end do
         151             :    end do
         152             :
         153      746136 :    call kiss_gen%finalize()
         154             :
         155      746136 : end subroutine mcica_subcol_lw
         156             :
         157             : !========================================================================================
         158             :
         159      389263 : subroutine mcica_subcol_sw( &
         160             :    kdist, nbnd, ngpt, ncol, nlay, nver, changeseed, &
         161      389263 :    pmid, cldfrac, tauc, ssac, asmc,     &
         162      389263 :    taucmcl, ssacmcl, asmcmcl)
         163             :
         164             :    ! Arrays use CAM vertical index convention: index increases from top to bottom.
         165             :    ! This index ordering is assumed in the maximum-random overlap algorithm which starts
         166             :    ! at the top of a column and marches down, with each layer depending on the state
         167             :    ! of the layer above it.
         168             :    !
         169             :    ! For GCM mode, changeseed must be offset between LW and SW by at least the
         170             :    ! number of subcolumns
         171             :
         172             :    ! arguments
         173             :    class(ty_gas_optics_rrtmgp), intent(in) :: kdist  ! spectral information
         174             :    integer,  intent(in)  :: nbnd                    ! number of spectral bands
         175             :    integer,  intent(in)  :: ngpt                    ! number of subcolumns (g-point intervals)
         176             :    integer,  intent(in)  :: ncol                    ! number of columns
         177             :    integer,  intent(in)  :: nlay                    ! number of vertical layers in radiation calc;
         178             :                                                     !   may include an "extra layer"
         179             :    integer,  intent(in)  :: nver                    ! number of CAM's vertical layers in rad calc
         180             :    integer,  intent(in)  :: changeseed              ! if the subcolumn generator is called multiple times,
         181             :                                                     ! permute the seed between each call.
         182             :    real(r8), intent(in)  :: pmid(ncol,nlay)         ! layer midpoint pressures (Pa)
         183             :    real(r8), intent(in)  :: cldfrac(ncol,nver)      ! layer cloud fraction
         184             :    real(r8), intent(in)  :: tauc(nbnd,ncol,nver)    ! cloud optical depth
         185             :    real(r8), intent(in)  :: ssac(nbnd,ncol,nver)    ! cloud single scattering albedo (non-delta scaled)
         186             :    real(r8), intent(in)  :: asmc(nbnd,ncol,nver)    ! cloud asymmetry parameter (non-delta scaled)
         187             :
         188             :
         189             :    real(r8), intent(out) :: taucmcl(ngpt,ncol,nver) ! subcolumn cloud optical depth [mcica]
         190             :    real(r8), intent(out) :: ssacmcl(ngpt,ncol,nver) ! subcolumn cloud single scattering albedo [mcica]
         191             :    real(r8), intent(out) :: asmcmcl(ngpt,ncol,nver) ! subcolumn cloud asymmetry parameter [mcica]
         192             :
         193             :    ! Local vars
         194             :
         195             :    integer :: i, isubcol, k, n
         196             :
         197             :    real(r8), parameter :: cldmin = 1.0e-80_r8  ! min cloud fraction
         198      778526 :    real(r8) :: cldf(ncol,nver)      ! cloud fraction clipped to cldmin
         199             :
         200      389263 :    type(ShrKissRandGen) :: kiss_gen  ! KISS RNG object
         201      778526 :    integer  :: kiss_seed(ncol,4)
         202      778526 :    real(r8) :: rand_num_1d(ncol,1)   ! random number (kissvec)
         203      778526 :    real(r8) :: rand_num(ncol,nver)   ! random number (kissvec)
         204             :
         205      778526 :    real(r8) :: cdf(ngpt,ncol,nver)   ! random numbers
         206      389263 :    logical  :: iscloudy(ngpt,ncol,nver)   ! flag that says whether a gridbox is cloudy
         207             :    !------------------------------------------------------------------------------------------
         208             :
         209             :    ! clip cloud fraction
         210   581226622 :    cldf(:,:) = cldfrac(:ncol,:)
         211   581226622 :    where (cldf(:,:) < cldmin)
         212             :       cldf(:,:) = 0._r8
         213             :    end where
         214             :
         215             :    ! Create a seed that depends on the state of the columns.
         216             :    ! Use pmid from bottom four layers.
         217     6245563 :    do i = 1, ncol
         218     5856300 :       kiss_seed(i,1) = (pmid(i,nlay)   - int(pmid(i,nlay)))    * 1000000000
         219     5856300 :       kiss_seed(i,2) = (pmid(i,nlay-1) - int(pmid(i,nlay-1)))  * 1000000000
         220     5856300 :       kiss_seed(i,3) = (pmid(i,nlay-2) - int(pmid(i,nlay-2)))  * 1000000000
         221     6245563 :       kiss_seed(i,4) = (pmid(i,nlay-3) - int(pmid(i,nlay-3)))  * 1000000000
         222             :    end do
         223             :
         224             :    ! create the RNG object
         225      389263 :    kiss_gen = ShrKissRandGen(kiss_seed)
         226             :
         227             :    ! Advance randum number generator by changeseed values
         228      778526 :    do i = 1, changeSeed
         229      778526 :       call kiss_gen%random(rand_num_1d)
         230             :    end do
         231             :
         232             :    ! Generate random numbers in each subcolumn at every level
         233    43986719 :    do isubcol = 1,ngpt
         234    43597456 :       call kiss_gen%random(rand_num)
         235 65097770927 :       cdf(isubcol,:,:) = rand_num(:,:)
         236             :    enddo
         237             :
         238             :    ! Maximum-Random overlap
         239             :    ! i) pick a random number for top layer.
         240             :    ! ii) walk down the column:
         241             :    !    - if the layer above is cloudy, use the same random number as in the layer above
         242             :    !    - if the layer above is clear, use a new random number
         243             :
         244    36201459 :    do k = 2, nver
         245   574981059 :       do i = 1, ncol
         246 60917906996 :          do isubcol = 1, ngpt
         247 60882094800 :             if (cdf(isubcol,i,k-1) > 1._r8 - cldf(i,k-1) ) then
         248  4754772523 :                cdf(isubcol,i,k) = cdf(isubcol,i,k-1)
         249             :             else
         250 55588542677 :                cdf(isubcol,i,k) = cdf(isubcol,i,k) * (1._r8 - cldf(i,k-1))
         251             :             end if
         252             :          end do
         253             :       end do
         254             :    end do
         255             :
         256    36590722 :    do k = 1, nver
         257 61580447422 :       iscloudy(:,:,k) = (cdf(:,:,k) >= 1._r8 - spread(cldf(:,k), dim=1, nCopies=ngpt) )
         258             :    end do
         259             :
         260             :    ! -- generate subcolumns for homogeneous clouds -----
         261             :    ! where there is a cloud, set the subcolumn cloud properties;
         262             :    ! incoming tauc should be in-cloud quantites and not grid-averaged quantities
         263    36590722 :    do k = 1,nver
         264   581226622 :       do i = 1,ncol
         265 61580058159 :          do isubcol = 1,ngpt
         266 61543856700 :             if (iscloudy(isubcol,i,k) .and. (cldf(i,k) > 0._r8) ) then
         267  4805438153 :                n = kdist%convert_gpt2band(isubcol)
         268  4805438153 :                taucmcl(isubcol,i,k) = tauc(n,i,k)
         269  4805438153 :                ssacmcl(isubcol,i,k) = ssac(n,i,k)
         270  4805438153 :                asmcmcl(isubcol,i,k) = asmc(n,i,k)
         271             :             else
         272 56193782647 :                taucmcl(isubcol,i,k) = 0._r8
         273 56193782647 :                ssacmcl(isubcol,i,k) = 1._r8
         274 56193782647 :                asmcmcl(isubcol,i,k) = 0._r8
         275             :             end if
         276             :          end do
         277             :       end do
         278             :    end do
         279             :
         280      389263 :    call kiss_gen%finalize()
         281             :
         282      389263 : end subroutine mcica_subcol_sw
         283             :
         284             :
         285             : end module mcica_subcol_gen
         286             :

Generated by: LCOV version 1.14