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