diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index efea61226..9e9f8c70c 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -956,7 +956,6 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & call progcld5 (plyr,plvl,tlyr,tvly,qlyr,qstl,rhly,tracer1, & ! --- inputs xlat,xlon,slmsk,dz,delp, & ntrac-1, ntcw-1,ntiw-1,ntrw-1, & -!mz ntsw-1,ntgl-1, & im, lmk, lmp, icloud, uni_cld, lmfshal, lmfdeep2, & cldcov(:,1:LMK),effrl_inout(:,:), & effri_inout(:,:), effrs_inout(:,:), & diff --git a/physics/HWRF_mcica_random_numbers.F90 b/physics/HWRF_mcica_random_numbers.F90 deleted file mode 100644 index b2f2d20dd..000000000 --- a/physics/HWRF_mcica_random_numbers.F90 +++ /dev/null @@ -1,109 +0,0 @@ - module mcica_random_numbers - - ! Generic module to wrap random number generators. - ! The module defines a type that identifies the particular stream of random - ! numbers, and has procedures for initializing it and getting real numbers - ! in the range 0 to 1. - ! This version uses the Mersenne Twister to generate random numbers on [0, 1]. - ! - use MersenneTwister, only: randomNumberSequence, & ! The random number engine. - new_RandomNumberSequence, getRandomReal -!! mji -!! use time_manager_mod, only: time_type, get_date - -!mz use parkind, only : im => kind_im, rb => kind_rb - use machine, only: im => kind_io4, rb => kind_phys - - implicit none - private - - type randomNumberStream - type(randomNumberSequence) :: theNumbers - end type randomNumberStream - - interface getRandomNumbers - module procedure getRandomNumber_Scalar, getRandomNumber_1D, getRandomNumber_2D - end interface getRandomNumbers - - interface initializeRandomNumberStream - module procedure initializeRandomNumberStream_S, initializeRandomNumberStream_V - end interface initializeRandomNumberStream - - public :: randomNumberStream, & - initializeRandomNumberStream, getRandomNumbers -!! mji -!! initializeRandomNumberStream, getRandomNumbers, & -!! constructSeed -contains - ! --------------------------------------------------------- - ! Initialization - ! --------------------------------------------------------- - function initializeRandomNumberStream_S(seed) result(new) - integer(kind=im), intent( in) :: seed - type(randomNumberStream) :: new - - new%theNumbers = new_RandomNumberSequence(seed) - - end function initializeRandomNumberStream_S - ! --------------------------------------------------------- - function initializeRandomNumberStream_V(seed) result(new) - integer(kind=im), dimension(:), intent( in) :: seed - type(randomNumberStream) :: new - - new%theNumbers = new_RandomNumberSequence(seed) - - end function initializeRandomNumberStream_V - - ! --------------------------------------------------------- - ! Procedures for drawing random numbers - ! --------------------------------------------------------- - subroutine getRandomNumber_Scalar(stream, number) - type(randomNumberStream), intent(inout) :: stream - real(kind=rb), intent( out) :: number - - number = getRandomReal(stream%theNumbers) - end subroutine getRandomNumber_Scalar - ! --------------------------------------------------------- - subroutine getRandomNumber_1D(stream, numbers) - type(randomNumberStream), intent(inout) :: stream - real(kind=rb), dimension(:), intent( out) :: numbers - - ! Local variables - integer(kind=im) :: i - - do i = 1, size(numbers) - numbers(i) = getRandomReal(stream%theNumbers) - end do - end subroutine getRandomNumber_1D - ! --------------------------------------------------------- - subroutine getRandomNumber_2D(stream, numbers) - type(randomNumberStream), intent(inout) :: stream - real(kind=rb), dimension(:, :), intent( out) :: numbers - - ! Local variables - integer(kind=im) :: i - - do i = 1, size(numbers, 2) - call getRandomNumber_1D(stream, numbers(:, i)) - end do - end subroutine getRandomNumber_2D - -! mji -! ! --------------------------------------------------------- -! ! Constructing a unique seed from grid cell index and model date/time -! ! Once we have the GFDL stuff we'll add the year, month, day, hour, minute -! ! --------------------------------------------------------- -! function constructSeed(i, j, time) result(seed) -! integer(kind=im), intent( in) :: i, j -! type(time_type), intent( in) :: time -! integer(kind=im), dimension(8) :: seed -! -! ! Local variables -! integer(kind=im) :: year, month, day, hour, minute, second -! -! -! call get_date(time, year, month, day, hour, minute, second) -! seed = (/ i, j, year, month, day, hour, minute, second /) -! end function constructSeed - - end module mcica_random_numbers diff --git a/physics/HWRF_mersenne_twister.F90 b/physics/HWRF_mersenne_twister.F90 deleted file mode 100644 index f9e3b0b0a..000000000 --- a/physics/HWRF_mersenne_twister.F90 +++ /dev/null @@ -1,304 +0,0 @@ -! Fortran-95 implementation of the Mersenne Twister 19937, following -! the C implementation described below (code mt19937ar-cok.c, dated 2002/2/10), -! adapted cosmetically by making the names more general. -! Users must declare one or more variables of type randomNumberSequence in the calling -! procedure which are then initialized using a required seed. If the -! variable is not initialized the random numbers will all be 0. -! For example: -! program testRandoms -! use RandomNumbers -! type(randomNumberSequence) :: randomNumbers -! integer :: i -! -! randomNumbers = new_RandomNumberSequence(seed = 100) -! do i = 1, 10 -! print ('(f12.10, 2x)'), getRandomReal(randomNumbers) -! end do -! end program testRandoms -! -! Fortran-95 implementation by -! Robert Pincus -! NOAA-CIRES Climate Diagnostics Center -! Boulder, CO 80305 -! email: Robert.Pincus@colorado.edu -! -! This documentation in the original C program reads: -! ------------------------------------------------------------- -! A C-program for MT19937, with initialization improved 2002/2/10. -! Coded by Takuji Nishimura and Makoto Matsumoto. -! This is a faster version by taking Shawn Cokus's optimization, -! Matthe Bellew's simplification, Isaku Wada's real version. -! -! Before using, initialize the state by using init_genrand(seed) -! or init_by_array(init_key, key_length). -! -! Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, -! All rights reserved. -! -! Redistribution and use in source and binary forms, with or without -! modification, are permitted provided that the following conditions -! are met: -! -! 1. Redistributions of source code must retain the above copyright -! notice, this list of conditions and the following disclaimer. -! -! 2. Redistributions in binary form must reproduce the above copyright -! notice, this list of conditions and the following disclaimer in the -! documentation and/or other materials provided with the distribution. -! -! 3. The names of its contributors may not be used to endorse or promote -! products derived from this software without specific prior written -! permission. -! -! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR -! A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR -! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -! EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR -! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF -! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING -! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -! -! -! Any feedback is very welcome. -! http://www.math.keio.ac.jp/matumoto/emt.html -! email: matumoto@math.keio.ac.jp -! ------------------------------------------------------------- - - module MersenneTwister -! ------------------------------------------------------------- - -!mz use parkind, only : im => kind_im, rb => kind_rb - use machine, only: im => kind_io4, rb => kind_phys - - implicit none - private - - ! Algorithm parameters - ! ------- - ! Period parameters - integer(kind=im), parameter :: blockSize = 624, & - M = 397, & - MATRIX_A = -1727483681, & ! constant vector a (0x9908b0dfUL) - UMASK = -2147483647-1, & ! most significant w-r bits (0x80000000UL) - LMASK = 2147483647 ! least significant r bits (0x7fffffffUL) - ! Tempering parameters - integer(kind=im), parameter :: TMASKB= -1658038656, & ! (0x9d2c5680UL) - TMASKC= -272236544 ! (0xefc60000UL) - ! ------- - - ! The type containing the state variable - type randomNumberSequence - integer(kind=im) :: currentElement ! = blockSize - integer(kind=im), dimension(0:blockSize -1) :: state ! = 0 - end type randomNumberSequence - - interface new_RandomNumberSequence - module procedure initialize_scalar, initialize_vector - end interface new_RandomNumberSequence - - - public :: randomNumberSequence - public :: new_RandomNumberSequence, finalize_RandomNumberSequence, & - getRandomInt, getRandomPositiveInt, getRandomReal -! ------------------------------------------------------------- -contains - ! ------------------------------------------------------------- - ! Private functions - ! --------------------------- - function mixbits(u, v) - integer(kind=im), intent( in) :: u, v - integer(kind=im) :: mixbits - - mixbits = ior(iand(u, UMASK), iand(v, LMASK)) - end function mixbits - ! --------------------------- - function twist(u, v) - integer(kind=im), intent( in) :: u, v - integer(kind=im) :: twist - - ! Local variable - integer(kind=im), parameter, dimension(0:1) :: t_matrix = (/ 0_im, MATRIX_A /) - - twist = ieor(ishft(mixbits(u, v), -1_im), t_matrix(iand(v, 1_im))) - twist = ieor(ishft(mixbits(u, v), -1_im), t_matrix(iand(v, 1_im))) - end function twist - ! --------------------------- - subroutine nextState(twister) - type(randomNumberSequence), intent(inout) :: twister - - ! Local variables - integer(kind=im) :: k - - do k = 0, blockSize - M - 1 - twister%state(k) = ieor(twister%state(k + M), & - twist(twister%state(k), twister%state(k + 1_im))) - end do - do k = blockSize - M, blockSize - 2 - twister%state(k) = ieor(twister%state(k + M - blockSize), & - twist(twister%state(k), twister%state(k + 1_im))) - end do - twister%state(blockSize - 1_im) = ieor(twister%state(M - 1_im), & - twist(twister%state(blockSize - 1_im), twister%state(0_im))) - twister%currentElement = 0_im - - end subroutine nextState - ! --------------------------- - elemental function temper(y) - integer(kind=im), intent(in) :: y - integer(kind=im) :: temper - - integer(kind=im) :: x - - ! Tempering - x = ieor(y, ishft(y, -11)) - x = ieor(x, iand(ishft(x, 7), TMASKB)) - x = ieor(x, iand(ishft(x, 15), TMASKC)) - temper = ieor(x, ishft(x, -18)) - end function temper - ! ------------------------------------------------------------- - ! Public (but hidden) functions - ! -------------------- - function initialize_scalar(seed) result(twister) - integer(kind=im), intent(in ) :: seed - type(randomNumberSequence) :: twister - - integer(kind=im) :: i - ! See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. In the previous versions, - ! MSBs of the seed affect only MSBs of the array state[]. - ! 2002/01/09 modified by Makoto Matsumoto - - twister%state(0) = iand(seed, -1_im) - do i = 1, blockSize - 1 ! ubound(twister%state) - twister%state(i) = 1812433253_im * ieor(twister%state(i-1), & - ishft(twister%state(i-1), -30_im)) + i - twister%state(i) = iand(twister%state(i), -1_im) ! for >32 bit machines - end do - twister%currentElement = blockSize - end function initialize_scalar - ! ------------------------------------------------------------- - function initialize_vector(seed) result(twister) - integer(kind=im), dimension(0:), intent(in) :: seed - type(randomNumberSequence) :: twister - - integer(kind=im) :: i, j, k, nFirstLoop, nWraps - - nWraps = 0 - twister = initialize_scalar(19650218_im) - - nFirstLoop = max(blockSize, size(seed)) - do k = 1, nFirstLoop - i = mod(k + nWraps, blockSize) - j = mod(k - 1, size(seed)) - if(i == 0) then - twister%state(i) = twister%state(blockSize - 1) - twister%state(1) = ieor(twister%state(1), & - ieor(twister%state(1-1), & - ishft(twister%state(1-1), -30_im)) * 1664525_im) + & - seed(j) + j ! Non-linear - twister%state(i) = iand(twister%state(i), -1_im) ! for >32 bit machines - nWraps = nWraps + 1 - else - twister%state(i) = ieor(twister%state(i), & - ieor(twister%state(i-1), & - ishft(twister%state(i-1), -30_im)) * 1664525_im) + & - seed(j) + j ! Non-linear - twister%state(i) = iand(twister%state(i), -1_im) ! for >32 bit machines - end if - end do - - ! - ! Walk through the state array, beginning where we left off in the block above - ! - do i = mod(nFirstLoop, blockSize) + nWraps + 1, blockSize - 1 - twister%state(i) = ieor(twister%state(i), & - ieor(twister%state(i-1), & - ishft(twister%state(i-1), -30_im)) * 1566083941_im) - i ! Non-linear - twister%state(i) = iand(twister%state(i), -1_im) ! for >32 bit machines - end do - - twister%state(0) = twister%state(blockSize - 1) - - do i = 1, mod(nFirstLoop, blockSize) + nWraps - twister%state(i) = ieor(twister%state(i), & - ieor(twister%state(i-1), & - ishft(twister%state(i-1), -30_im)) * 1566083941_im) - i ! Non-linear - twister%state(i) = iand(twister%state(i), -1_im) ! for >32 bit machines - end do - - twister%state(0) = UMASK - twister%currentElement = blockSize - - end function initialize_vector - ! ------------------------------------------------------------- - ! Public functions - ! -------------------- - function getRandomInt(twister) - type(randomNumberSequence), intent(inout) :: twister - integer(kind=im) :: getRandomInt - ! Generate a random integer on the interval [0,0xffffffff] - ! Equivalent to genrand_int32 in the C code. - ! Fortran doesn't have a type that's unsigned like C does, - ! so this is integers in the range -2**31 - 2**31 - ! All functions for getting random numbers call this one, - ! then manipulate the result - - if(twister%currentElement >= blockSize) call nextState(twister) - - getRandomInt = temper(twister%state(twister%currentElement)) - twister%currentElement = twister%currentElement + 1 - - end function getRandomInt - ! -------------------- - function getRandomPositiveInt(twister) - type(randomNumberSequence), intent(inout) :: twister - integer(kind=im) :: getRandomPositiveInt - ! Generate a random integer on the interval [0,0x7fffffff] - ! or [0,2**31] - ! Equivalent to genrand_int31 in the C code. - - ! Local integers - integer(kind=im) :: localInt - - localInt = getRandomInt(twister) - getRandomPositiveInt = ishft(localInt, -1) - - end function getRandomPositiveInt - ! -------------------- - ! -------------------- -!! mji - modified Jan 2007, double converted to rrtmg real kind type - function getRandomReal(twister) - type(randomNumberSequence), intent(inout) :: twister -! double precision :: getRandomReal - real(kind=rb) :: getRandomReal - ! Generate a random number on [0,1] - ! Equivalent to genrand_real1 in the C code - ! The result is stored as double precision but has 32 bit resolution - - integer(kind=im) :: localInt - - localInt = getRandomInt(twister) - if(localInt < 0) then -! getRandomReal = dble(localInt + 2.0d0**32)/(2.0d0**32 - 1.0d0) - getRandomReal = (localInt + 2.0**32_rb)/(2.0**32_rb - 1.0_rb) - else -! getRandomReal = dble(localInt )/(2.0d0**32 - 1.0d0) - getRandomReal = (localInt )/(2.0**32_rb - 1.0_rb) - end if - - end function getRandomReal - ! -------------------- - subroutine finalize_RandomNumberSequence(twister) - type(randomNumberSequence), intent(inout) :: twister - - twister%currentElement = blockSize - twister%state(:) = 0_im - end subroutine finalize_RandomNumberSequence - - ! -------------------- - - end module MersenneTwister - diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index 36c928b2d..119a2d1d0 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -2753,63 +2753,6 @@ subroutine progcld5 & enddo enddo endif -!mz - if (icloud .ne. 0) then -! assign/calculate efective radii for cloud water, ice, rain, snow - - do k = 1, NLAY - do i = 1, IX - rew(i,k) = reliq_def ! default liq radius to 10 micron - rei(i,k) = reice_def ! default ice radius to 50 micron - rer(i,k) = rrain_def ! default rain radius to 1000 micron - res(i,k) = rsnow_def ! default snow radius to 250 micron - enddo - enddo -!> -# Compute effective liquid cloud droplet radius over land. - do i = 1, IX - if (nint(slmsk(i)) == 1) then - do k = 1, NLAY - tem1 = min(1.0, max(0.0, (con_ttp-tlyr(i,k))*0.05)) - rew(i,k) = 5.0 + 5.0 * tem1 - enddo - endif - enddo - -!> -# Compute effective ice cloud droplet radius following Heymsfield -!! and McFarquhar (1996) \cite heymsfield_and_mcfarquhar_1996. - - do k = 1, NLAY - do i = 1, IX - tem2 = tlyr(i,k) - con_ttp - - if (cip(i,k) > 0.0) then - tem3 = gord * cip(i,k) * plyr(i,k) / (delp(i,k)*tvly(i,k)) - - if (tem2 < -50.0) then - rei(i,k) = (1250.0/9.917) * tem3 ** 0.109 - elseif (tem2 < -40.0) then - rei(i,k) = (1250.0/9.337) * tem3 ** 0.08 - elseif (tem2 < -30.0) then - rei(i,k) = (1250.0/9.208) * tem3 ** 0.055 - else - rei(i,k) = (1250.0/9.387) * tem3 ** 0.031 - endif - rei(i,k) = max(25.,rei(i,k)) !mz* HWRF - endif - rei(i,k) = min(rei(i,k), 135.72) !- 1.0315*rei<= 140 microns - enddo - enddo - -!mz -!> -# Compute effective snow cloud droplet radius - do k = 1, NLAY - do i = 1, IX - res(i,k) = 10.0 - enddo - enddo - - endif ! end icloud -!mz end do k = 1, NLAY do i = 1, IX clouds(i,k,1) = cldtot(i,k) diff --git a/physics/radlw_main.F90 b/physics/radlw_main.F90 index 89609c283..14a28cf6b 100644 --- a/physics/radlw_main.F90 +++ b/physics/radlw_main.F90 @@ -284,7 +284,6 @@ module rrtmg_lw & con_amw, con_amo3 use mersenne_twister, only : random_setseed, random_number, & & random_stat -!mz use machine, only : kind_phys, & & im => kind_io4, rb => kind_phys @@ -635,31 +634,6 @@ subroutine rrtmg_lw_run & real (kind=kind_phys), dimension(:,:,:),intent(in):: & & aeraod, aerssa -!mz* HWRF -- OUTPUT from mcica_subcol_lw - real(kind=kind_phys),dimension(ngptlw,npts,nlay) :: cldfmcl ! Cloud fraction - ! Dimensions: (ngptlw,ncol,nlay) - real(kind=kind_phys),dimension(ngptlw,npts,nlay) :: ciwpmcl ! In-cloud ice water path (g/m2) - ! Dimensions: (ngptlw,ncol,nlay) - real(kind=kind_phys),dimension(ngptlw,npts,nlay) :: clwpmcl ! In-cloud liquid water path (g/m2) - ! Dimensions: (ngptlw,ncol,nlay) - real(kind=kind_phys),dimension(ngptlw,npts,nlay) :: cswpmcl ! In-cloud snow water path (g/m2) - ! Dimensions: (ngptlw,ncol,nlay) - real(kind=kind_phys),dimension(npts,nlay) :: relqmcl ! Cloud water drop effective radius (microns) - ! Dimensions: (ncol,nlay) - real(kind=kind_phys),dimension(npts,nlay) :: reicmcl ! Cloud ice effective size (microns) - ! Dimensions: (ncol,nlay) - real(kind=kind_phys),dimension(npts,nlay) :: resnmcl ! Snow effective size (microns) - ! Dimensions: (ncol,nlay) - real(kind=kind_phys),dimension(ngptlw,npts,nlay) :: taucmcl ! In-cloud optical depth - ! Dimensions: (ngptlw,ncol,nlay) -! real(kind=kind_phys),dimension(npts,nlay,nbands) :: tauaer ! Aerosol optical depth -! ! Dimensions: (ncol,nlay,nbndlw) -!mz* output from cldprmc - integer :: ncbands ! number of cloud spectral bands - real(kind=kind_phys),dimension(ngptlw,nlay) :: taucmc ! cloud optical depth [mcica] - ! Dimensions: (ngptlw,nlayers) -!mz - ! --- outputs: real (kind=kind_phys), dimension(:,:), intent(inout) :: hlwc real (kind=kind_phys), dimension(:,:), intent(inout) :: & @@ -681,11 +655,6 @@ subroutine rrtmg_lw_run & logical, intent(in) :: lslwr ! --- locals: -! mz* - Add height of each layer for exponential-random cloud overlap -! This will be derived below from the dzlyr in each layer - real (kind=kind_phys), dimension( npts,nlay ) :: hgt - real (kind=kind_phys) :: dzsum - real (kind=kind_phys), dimension(0:nlp1) :: cldfrc real (kind=kind_phys), dimension(0:nlay) :: totuflux, totdflux, & @@ -699,11 +668,6 @@ subroutine rrtmg_lw_run & & selffac, selffrac, forfac, forfrac, minorfrac, scaleminor, & & scaleminorn2, temcol, dz -!mz* - real(kind=rb),dimension(0:nlay,nbands) :: planklay,planklev - real(kind=rb),dimension(0:nlay) :: pz - -! real(kind=rb) :: plankbnd(nbndlw) real (kind=kind_phys), dimension(nbands,0:nlay) :: pklev, pklay real (kind=kind_phys), dimension(nlay,nbands) :: htrb @@ -711,26 +675,7 @@ subroutine rrtmg_lw_run & real (kind=kind_phys), dimension(nbands,npts,nlay) :: taucld3 real (kind=kind_phys), dimension(ngptlw,nlay) :: fracs, tautot real (kind=kind_phys), dimension(nlay,ngptlw) :: fracs_r -!mz rtrnmc_mcica - real (kind=kind_phys), dimension(nlay,ngptlw) :: taut -!mz* Atmosphere/clouds - cldprop - real(kind=kind_phys), dimension(ngptlw,nlay) :: cldfmc, & - & cldfmc_save ! cloud fraction [mcica] - ! Dimensions: (ngptlw,nlay) - real(kind=kind_phys), dimension(ngptlw,nlay) :: ciwpmc ! in-cloud ice water path [mcica] - ! Dimensions: (ngptlw,nlay) - real(kind=kind_phys), dimension(ngptlw,nlay) :: clwpmc ! in-cloud liquid water path [mcica] - ! Dimensions: (ngptlw,nlay) - real(kind=kind_phys), dimension(ngptlw,nlay) :: cswpmc ! in-cloud snow path [mcica] - ! Dimensions: (ngptlw,nlay) - real(kind=kind_phys), dimension(nlay) :: relqmc ! liquid particle effective radius (microns) - ! Dimensions: (nlay) - real(kind=kind_phys), dimension(nlay) :: reicmc ! ice particle effective size (microns) - ! Dimensions: (nlay) - real(kind=kind_phys), dimension(nlay) :: resnmc ! snow effective size (microns) - ! Dimensions: (nlay) - - + real(kind=kind_phys), dimension(ngptlw,nlay) :: cldfmc real (kind=kind_phys), dimension(nbands) :: semiss, secdiff ! --- column amount of absorbing gases: @@ -752,8 +697,7 @@ subroutine rrtmg_lw_run & integer, dimension(npts) :: ipseed integer, dimension(nlay) :: jp, jt, jt1, indself, indfor, indminor integer :: laytrop, iplon, i, j, k, k1 - ! mz* added local arrays for RRTMG - integer :: irng, permuteseed,ig + integer :: ig integer :: inflglw, iceflglw, liqflglw logical :: lcf1 integer :: istart ! beginning band of calculation @@ -850,43 +794,6 @@ subroutine rrtmg_lw_run & stemp = sfgtmp(iplon) ! surface ground temp if (iovr == 3) delgth= de_lgth(iplon) ! clouds decorr-length -! mz*: HWRF - if (iovr == 4 ) then - -!Add layer height needed for exponential (icld=4) and -! exponential-random (icld=5) overlap options - - !iplon = 1 - irng = 0 - permuteseed = 150 - -!mz* Derive height - dzsum =0.0 - do k = 1,nlay - hgt(iplon,k)= dzsum+0.5*dzlyr(iplon,k)*1000. !km->m - dzsum = dzsum+ dzlyr(iplon,k)*1000. - enddo - -! Zero out cloud optical properties here; not used when passing physical properties -! to radiation and taucld is calculated in radiation - do k = 1, nlay - do j = 1, nbands - taucld3(j,iplon,k) = 0.0 - enddo - enddo - - call mcica_subcol_lw(1, iplon, nlay, iovr, permuteseed, & - & irng, plyr, hgt, & - & cld_cf, cld_iwp, cld_lwp,cld_swp, & - & cld_ref_ice, cld_ref_liq, & - & cld_ref_snow, taucld3, & - & cldfmcl, & !--output - & ciwpmcl, clwpmcl, cswpmcl, reicmcl, relqmcl, & - & resnmcl, taucmcl) - - endif -!mz* end - !> -# Prepare atmospheric profile for use in rrtm. ! the vertical index of internal array is from surface to top @@ -987,23 +894,6 @@ subroutine rrtmg_lw_run & cda3(k) = cld_swp(iplon,k1) cda4(k) = cld_ref_snow(iplon,k1) enddo - ! HWRF RRMTG - if (iovr == 4) then !mz HWRF - do k = 1, nlay - k1 = nlp1 - k - do ig = 1, ngptlw - cldfmc(ig,k) = cldfmcl(ig,iplon,k1) - taucmc(ig,k) = taucmcl(ig,iplon,k1) - ciwpmc(ig,k) = ciwpmcl(ig,iplon,k1) - clwpmc(ig,k) = clwpmcl(ig,iplon,k1) - !mz cswpmc(ig,k) = cswpmcl(ig,iplon,k1) - cswpmc(ig,k) = 0.0 - enddo - reicmc(k) = reicmcl(iplon,k1) - relqmc(k) = relqmcl(iplon,k1) - resnmc(k) = resnmcl(iplon,k1) - enddo - endif else ! use diagnostic cloud method do k = 1, nlay k1 = nlp1 - k @@ -1111,24 +1001,6 @@ subroutine rrtmg_lw_run & cda3(k) = cld_swp(iplon,k) cda4(k) = cld_ref_snow(iplon,k) enddo - if (iovr == 4) then -!mz* Move incoming GCM cloud arrays to RRTMG cloud arrays. -!For GCM input, incoming reicmcl is defined based on selected -!ice parameterization (inflglw) - do k = 1, nlay - do ig = 1, ngptlw - cldfmc(ig,k) = cldfmcl(ig,iplon,k) - taucmc(ig,k) = taucmcl(ig,iplon,k) - ciwpmc(ig,k) = ciwpmcl(ig,iplon,k) - clwpmc(ig,k) = clwpmcl(ig,iplon,k) - !mz cswpmc(ig,k) = cswpmcl(ig,iplon,k) - cswpmc(ig,k) = 0.0 - enddo - reicmc(k) = reicmcl(iplon,k) - relqmc(k) = relqmcl(iplon,k) - resnmc(k) = resnmcl(iplon,k) - enddo - endif else ! use diagnostic cloud method do k = 1, nlay cldfrc(k)= cld_cf(iplon,k) @@ -1204,15 +1076,6 @@ subroutine rrtmg_lw_run & if ( lcf1 ) then - !mz* for HWRF, save cldfmc with mcica - if (iovr == 4) then - do k = 1, nlay - do ig = 1, ngptlw - cldfmc_save(ig,k)=cldfmc (ig,k) - enddo - enddo - endif - call cldprop & ! --- inputs: & ( cldfrc,clwp,relw,ciwp,reiw,cda1,cda2,cda3,cda4, & @@ -1221,15 +1084,6 @@ subroutine rrtmg_lw_run & & cldfmc, taucld & & ) - if (iovr == 4) then - !mz for HWRF, still using mcica cldfmc - do k = 1, nlay - do ig = 1, ngptlw - cldfmc(ig,k)=cldfmc_save(ig,k) - enddo - enddo - endif - ! --- ... save computed layer cloud optical depth for output ! rrtm band-7 is apprx 10mu channel (or use spectral mean of bands 6-8) @@ -1249,16 +1103,6 @@ subroutine rrtmg_lw_run & taucld = f_zero endif -!mz* HWRF: calculate taucmc with mcica - if (iovr == 4) then - call cldprmc(nlay, inflglw, iceflglw, liqflglw, & - & cldfmc, ciwpmc, & - & clwpmc, cswpmc, reicmc, relqmc, resnmc, & - & ncbands, taucmc, errmsg, errflg) - ! return immediately if cldprmc throws an error - if (errflg/=0) return - endif - ! if (lprnt) then ! print *,' after cldprop' ! print *,' clwp',clwp @@ -1992,8 +1836,6 @@ subroutine cldprop & ! --- ... call sub-column cloud generator -!mz* - if (iovr .ne. 4) then call mcica_subcol & ! --- inputs: & ( cldf, nlay, ipseed, dz, de_lgth, alpha, & @@ -2010,7 +1852,6 @@ subroutine cldprop & endif enddo enddo - endif !iovr endif ! end if_isubclw_block @@ -2243,39 +2084,39 @@ subroutine mcica_subcol & ! --- setup 2 sets of random numbers - call random_number ( rand2d, stat ) +! call random_number ( rand2d, stat ) - k1 = 0 - do k = 1, nlay - do n = 1, ngptlw - k1 = k1 + 1 - cdfunc(n,k) = rand2d(k1) - enddo - enddo +! k1 = 0 +! do k = 1, nlay +! do n = 1, ngptlw +! k1 = k1 + 1 +! cdfunc(n,k) = rand2d(k1) +! enddo +! enddo - call random_number ( rand2d, stat ) +! call random_number ( rand2d, stat ) - k1 = 0 - do k = 1, nlay - do n = 1, ngptlw - k1 = k1 + 1 - cdfun2(n,k) = rand2d(k1) - enddo - enddo +! k1 = 0 +! do k = 1, nlay +! do n = 1, ngptlw +! k1 = k1 + 1 +! cdfun2(n,k) = rand2d(k1) +! enddo +! enddo ! --- then working upward from the surface: ! if a random number (from an independent set: cdfun2) is smaller than ! alpha, then use the previous layer's number, otherwise use a new random ! number (keep the originally assigned one in cdfunc for that layer). - do k = 2, nlay - k1 = k - 1 - do n = 1, ngptlw - if ( cdfun2(n,k) < alpha(k) ) then - cdfunc(n,k) = cdfunc(n,k1) - endif - enddo - enddo +! do k = 2, nlay +! k1 = k - 1 +! do n = 1, ngptlw +! if ( cdfun2(n,k) < alpha(k) ) then +! cdfunc(n,k) = cdfunc(n,k1) +! endif +! enddo +! enddo end select @@ -7040,1023 +6881,6 @@ end subroutine taugb16 ! .................................. end subroutine taumol !! @} -!----------------------------------- - -!mz* exponential cloud overlapping subroutines -!------------------------------------------------------------------ -! Public subroutines -!------------------------------------------------------------------ -! mz* - Add height needed for exponential and exponential-random cloud overlap methods (icld=4 and 5, respectively) - subroutine mcica_subcol_lw(iplon, ncol, nlay, icld, permuteseed, & - & irng, play, hgt, & - & cldfrac, ciwp, clwp, cswp, rei, rel, res, tauc, & - & cldfmcl, & - & ciwpmcl, clwpmcl, cswpmcl, reicmcl, relqmcl, & - & resnmcl, taucmcl) - - use machine, only : im => kind_io4, rb => kind_phys -! ----- Input ----- -! Control - integer(kind=im), intent(in) :: iplon ! column/longitude index - integer(kind=im), intent(in) :: ncol ! number of columns - integer(kind=im), intent(in) :: nlay ! number of model layers - integer(kind=im), intent(in) :: icld ! clear/cloud, cloud overlap flag - integer(kind=im), intent(in) :: permuteseed ! if the cloud generator is called multiple times, - ! permute the seed between each call. - ! between calls for LW and SW, recommended - ! permuteseed differes by 'ngpt' - integer(kind=im), intent(inout) :: irng ! flag for random number generator - ! 0 = kissvec - ! 1 = Mersenne - ! Twister - -! Atmosphere - real(kind=rb), intent(in) :: play(:,:) ! layer pressures (mb) - ! Dimensions: (ncol,nlay) - -! mji - Add height - real(kind=rb), intent(in) :: hgt(:,:) ! layer height (m) - ! Dimensions: (ncol,nlay) - -! Atmosphere/clouds - cldprop - real(kind=rb), intent(in) :: cldfrac(:,:) ! layer cloud fraction - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: tauc(:,:,:) ! in-cloud optical depth - ! Dimensions: (nbndlw,ncol,nlay) -! real(kind=rb), intent(in) :: ssac(:,:,:) ! in-cloud single scattering albedo - ! Dimensions: (nbndlw,ncol,nlay) -! real(kind=rb), intent(in) :: asmc(:,:,:) ! in-cloud asymmetry parameter - ! Dimensions: (nbndlw,ncol,nlay) - real(kind=rb), intent(in) :: ciwp(:,:) ! in-cloud ice water path - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: clwp(:,:) ! in-cloud liquid water path - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: cswp(:,:) ! in-cloud snow path - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: rei(:,:) ! cloud ice particle size - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: rel(:,:) ! cloud liquid particle size - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: res(:,:) ! snow particle size - ! Dimensions: (ncol,nlay) - -! ----- Output ----- -! Atmosphere/clouds - cldprmc [mcica] - real(kind=rb), intent(out) :: cldfmcl(:,:,:) ! cloud fraction [mcica] - ! Dimensions: (ngptlw,ncol,nlay) - real(kind=rb), intent(out) :: ciwpmcl(:,:,:) ! in-cloud ice water path [mcica] - ! Dimensions: (ngptlw,ncol,nlay) - real(kind=rb), intent(out) :: clwpmcl(:,:,:) ! in-cloud liquid water path [mcica] - ! Dimensions: (ngptlw,ncol,nlay) - real(kind=rb), intent(out) :: cswpmcl(:,:,:) ! in-cloud snow path [mcica] - ! Dimensions: (ngptlw,ncol,nlay) - real(kind=rb), intent(out) :: relqmcl(:,:) ! liquid particle size (microns) - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(out) :: reicmcl(:,:) ! ice partcle size (microns) - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(out) :: resnmcl(:,:) ! snow partcle size (microns) - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(out) :: taucmcl(:,:,:) ! in-cloud optical depth [mcica] -!mz* - ! Dimensions: (ngptlw,ncol,nlay) -! real(kind=rb), intent(out) :: ssacmcl(:,:,:) ! in-cloud single scattering albedo [mcica] - ! Dimensions: (ngptlw,ncol,nlay) -! real(kind=rb), intent(out) :: asmcmcl(:,:,:) ! in-cloud asymmetry parameter [mcica] - ! Dimensions: (ngptlw,ncol,nlay) -! ----- Local ----- - -! Stochastic cloud generator variables [mcica] - integer(kind=im), parameter :: nsubclw = ngptlw ! number of sub-columns (g-point intervals) - integer(kind=im) :: ilev ! loop index - - real(kind=rb) :: pmid(ncol, nlay) ! layer pressures (Pa) -! real(kind=rb) :: pdel(ncol, nlay) ! layer pressure thickness (Pa) -! real(kind=rb) :: qi(ncol, nlay) ! ice water (specific humidity) -! real(kind=rb) :: ql(ncol, nlay) ! liq water (specific humidity) - -! Return if clear sky - if (icld.eq.0) return - -! NOTE: For GCM mode, permuteseed must be offset between LW and SW by at least the number of subcolumns - - -! Pass particle sizes to new arrays, no subcolumns for these properties yet -! Convert pressures from mb to Pa - - reicmcl(:ncol,:nlay) = rei(:ncol,:nlay) - relqmcl(:ncol,:nlay) = rel(:ncol,:nlay) - resnmcl(:ncol,:nlay) = res(:ncol,:nlay) - pmid(:ncol,:nlay) = play(:ncol,:nlay)*1.e2_rb - -! Generate the stochastic subcolumns of cloud optical properties for -! the longwave - call generate_stochastic_clouds (ncol, nlay, nsubclw, icld, irng, & - & pmid, hgt, cldfrac, clwp, ciwp, cswp, tauc, & - & cldfmcl, clwpmcl, ciwpmcl, cswpmcl, & - & taucmcl, permuteseed) - - end subroutine mcica_subcol_lw -!------------------------------------------------------------------------------------------------- - subroutine generate_stochastic_clouds(ncol, nlay, nsubcol, icld, & - & irng, pmid, hgt, cld, clwp, ciwp, cswp, tauc, & - & cld_stoch, clwp_stoch, ciwp_stoch, & - & cswp_stoch, tauc_stoch, changeSeed) -!------------------------------------------------------------------------------------------------- -!------------------------------------------------------------------------------------------------- -! Contact: Cecile Hannay (hannay@ucar.edu) -! -! Original code: Based on Raisanen et al., QJRMS, 2004. -! -! Modifications: -! 1) Generalized for use with RRTMG and added Mersenne Twister as the default -! random number generator, which can be changed to the optional kissvec random number generator -! with flag 'irng'. Some extra functionality has been commented or removed. -! Michael J. Iacono, AER, Inc., February 2007 -! 2) Activated exponential and exponential/random cloud overlap method -! Michael J. Iacono, AER, November 2017 -! -! Given a profile of cloud fraction, cloud water and cloud ice, we produce a set of subcolumns. -! Each layer within each subcolumn is homogeneous, with cloud fraction equal to zero or one -! and uniform cloud liquid and cloud ice concentration. -! The ensemble as a whole reproduces the probability function of cloud liquid and ice within each layer -! and obeys an overlap assumption in the vertical. -! -! Overlap assumption: -! The cloud are consistent with 5 overlap assumptions: random, maximum, maximum-random, exponential and exponential random. -! The default option is maximum-random (option 2) -! The options are: 1=random overlap, 2=max/random, 3=maximum overlap, 4=exponential overlap, 5=exp/random -! This is set with the variable "overlap" -! The exponential overlap uses also a length scale, Zo. (real, parameter :: Zo = 2500. ) -! -! Seed: -! If the stochastic cloud generator is called several times during the same timestep, -! one should change the seed between the call to insure that the -! subcolumns are different. -! This is done by changing the argument 'changeSeed' -! For example, if one wants to create a set of columns for the -! shortwave and another set for the longwave , -! use 'changeSeed = 1' for the first call and'changeSeed = 2' for the second call - -! PDF assumption: -! We can use arbitrary complicated PDFS. -! In the present version, we produce homogeneuous clouds (the simplest case). -! Future developments include using the PDF scheme of Ben Johnson. -! -! History file: -! Option to add diagnostics variables in the history file. (using FINCL in the namelist) -! nsubcol = number of subcolumns -! overlap = overlap type (1-3) -! Zo = length scale -! CLOUD_S = mean of the subcolumn cloud fraction ('_S" means Stochastic) -! CLDLIQ_S = mean of the subcolumn cloud water -! CLDICE_S = mean of the subcolumn cloud ice -! -! Note: -! Here: we force that the cloud condensate to be consistent with the cloud fraction -! i.e we only have cloud condensate when the cell is cloudy. -! In CAM: The cloud condensate and the cloud fraction are obtained from 2 different equations -! and the 2 quantities can be inconsistent (i.e. CAM can produce cloud fraction -! without cloud condensate or the opposite). -!----------------------------------------------------------------- - - use mcica_random_numbers -! The Mersenne Twister random number engine - use MersenneTwister, only: randomNumberSequence, & - & new_RandomNumberSequence, getRandomReal - use machine ,only : im => kind_io4, rb => kind_phys - - type(randomNumberSequence) :: randomNumbers - -! -- Arguments - - integer(kind=im), intent(in) :: ncol ! number of columns - integer(kind=im), intent(in) :: nlay ! number of layers - integer(kind=im), intent(in) :: icld ! clear/cloud, cloud overlap flag - integer(kind=im), intent(inout) :: irng ! flag for random number generator - ! 0 = kissvec - ! 1 = Mersenne Twister - integer(kind=im), intent(in) :: nsubcol ! number of sub-columns (g-point intervals) - integer(kind=im), optional, intent(in) :: changeSeed ! allows permuting seed - -! Column state (cloud fraction, cloud water, cloud ice) + variables needed to read physics state - real(kind=rb), intent(in) :: pmid(:,:) ! layer pressure (Pa) - ! Dimensions: (ncol,nlay) - - real(kind=rb), intent(in) :: hgt(:,:) ! layer height (m) - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: cld(:,:) ! cloud fraction - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: clwp(:,:) ! in-cloud liquid water path - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: ciwp(:,:) ! in-cloud ice water path - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: cswp(:,:) ! in-cloud snow path - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: tauc(:,:,:) ! in-cloud optical depth - ! Dimensions:(nbndlw,ncol,nlay) -! real(kind=rb), intent(in) :: ssac(:,:,:) ! in-cloud single scattering albedo - ! Dimensions: (nbndlw,ncol,nlay) - ! inactive - for future expansion -! real(kind=rb), intent(in) :: asmc(:,:,:) ! in-cloud asymmetry parameter - ! Dimensions: (nbndlw,ncol,nlay) - ! inactive - for future expansion - - real(kind=rb), intent(out) :: cld_stoch(:,:,:) ! subcolumn cloud fraction - ! Dimensions: (ngptlw,ncol,nlay) - real(kind=rb), intent(out) :: clwp_stoch(:,:,:) ! subcolumn in-cloud liquid water path - ! Dimensions: (ngptlw,ncol,nlay) - real(kind=rb), intent(out) :: ciwp_stoch(:,:,:) ! subcolumn in-cloud ice water path - ! Dimensions: (ngptlw,ncol,nlay) - real(kind=rb), intent(out) :: cswp_stoch(:,:,:) ! subcolumn in-cloud snow path - ! Dimensions: (ngptlw,ncol,nlay) - real(kind=rb), intent(out) :: tauc_stoch(:,:,:) ! subcolumn in-cloud optical depth - ! Dimensions: (ngptlw,ncol,nlay) -! real(kind=rb), intent(out) :: ssac_stoch(:,:,:)! subcolumn in-cloud single scattering albedo - ! Dimensions: (ngptlw,ncol,nlay) - ! inactive - for future expansion -! real(kind=rb), intent(out) :: asmc_stoch(:,:,:)! subcolumn in-cloud asymmetry parameter - ! Dimensions: (ngptlw,ncol,nlay) - ! inactive - for future expansion - -! -- Local variables - real(kind=rb) :: cldf(ncol,nlay) ! cloud fraction - -! Mean over the subcolumns (cloud fraction, cloud water , cloud ice) - inactive -! real(kind=rb) :: mean_cld_stoch(ncol, nlay) ! cloud fraction -! real(kind=rb) :: mean_clwp_stoch(ncol, nlay) ! cloud water -! real(kind=rb) :: mean_ciwp_stoch(ncol, nlay) ! cloud ice -! real(kind=rb) :: mean_tauc_stoch(ncol, nlay) ! cloud optical depth -! real(kind=rb) :: mean_ssac_stoch(ncol, nlay) ! cloud single scattering albedo -! real(kind=rb) :: mean_asmc_stoch(ncol, nlay) ! cloud asymmetry parameter - -! Set overlap - integer(kind=im) :: overlap ! 1 = random overlap, 2 = maximum-random, - ! 3 = maximum overlap, 4 = exponential, - ! 5 = exponential-random - real(kind=rb), parameter :: Zo = 2500._rb ! length scale (m) - real(kind=rb), dimension(ncol,nlay) :: alpha ! overlap parameter - -! Constants (min value for cloud fraction and cloud water and ice) - real(kind=rb), parameter :: cldmin = 1.0e-20_rb ! min cloud fraction -! real(kind=rb), parameter :: qmin = 1.0e-10_rb ! min cloud water and cloud ice (not used) - -! Variables related to random number and seed - real(kind=rb), dimension(nsubcol, ncol, nlay) :: CDF, CDF2 !random numbers - integer(kind=im), dimension(ncol) :: seed1, seed2, seed3, seed4 !seed to create random number (kissvec) - real(kind=rb), dimension(ncol) :: rand_num ! random number (kissvec) - integer(kind=im) :: iseed ! seed to create random number (Mersenne Teister) - real(kind=rb) :: rand_num_mt ! random number (Mersenne Twister) - -! Flag to identify cloud fraction in subcolumns - logical, dimension(nsubcol, ncol, nlay) :: iscloudy ! flag that says whether a gridbox is cloudy - -! Indices - integer(kind=im) :: ilev, isubcol, i, n ! indices - -!------------------------------------------------------------------- - -! Check that irng is in bounds; if not, set to default - if (irng .ne. 0) irng = 1 - -! Pass input cloud overlap setting to local variable - overlap = icld - -! Ensure that cloud fractions are in bounds - do ilev = 1, nlay - do i = 1, ncol - cldf(i,ilev) = cld(i,ilev) - if (cldf(i,ilev) < cldmin) then - cldf(i,ilev) = 0._rb - endif - enddo - enddo - -! ----- Create seed -------- - -! Advance randum number generator by changeseed values - if (irng.eq.0) then -! For kissvec, create a seed that depends on the state of the columns. Maybe not the best way, but it works. -! Must use pmid from bottom four layers. - do i=1,ncol - if (pmid(i,1).lt.pmid(i,2)) then - stop 'MCICA_SUBCOL: KISSVEC SEED GENERATOR REQUIRES PMID & - & FROM BOTTOM FOUR LAYERS.' - endif - seed1(i) = (pmid(i,1) - int(pmid(i,1))) * 1000000000_im - seed2(i) = (pmid(i,2) - int(pmid(i,2))) * 1000000000_im - seed3(i) = (pmid(i,3) - int(pmid(i,3))) * 1000000000_im - seed4(i) = (pmid(i,4) - int(pmid(i,4))) * 1000000000_im - enddo - do i=1,changeSeed - call kissvec(seed1, seed2, seed3, seed4, rand_num) - enddo - elseif (irng.eq.1) then - randomNumbers = new_RandomNumberSequence(seed = changeSeed) - endif - -! ------ Apply overlap assumption -------- - -! generate the random numbers - - select case (overlap) - - case(1) -! Random overlap -! i) pick a random value at every level - - if (irng.eq.0) then - do isubcol = 1,nsubcol - do ilev = 1,nlay - call kissvec(seed1, seed2, seed3, seed4, rand_num) ! we get different random number for each level - CDF(isubcol,:,ilev) = rand_num - enddo - enddo - elseif (irng.eq.1) then - do isubcol = 1, nsubcol - do i = 1, ncol - do ilev = 1, nlay - rand_num_mt = getRandomReal(randomNumbers) - CDF(isubcol,i,ilev) = rand_num_mt - enddo - enddo - enddo - endif - - case(2) -! Maximum-Random overlap -! i) pick a random number for top layer. -! ii) walk down the column: -! - if the layer above is cloudy, we use the same random number than in the layer above -! - if the layer above is clear, we use a new random number - - if (irng.eq.0) then - do isubcol = 1,nsubcol - do ilev = 1,nlay - call kissvec(seed1, seed2, seed3, seed4, rand_num) - CDF(isubcol,:,ilev) = rand_num - enddo - enddo - elseif (irng.eq.1) then - do isubcol = 1, nsubcol - do i = 1, ncol - do ilev = 1, nlay - rand_num_mt = getRandomReal(randomNumbers) - CDF(isubcol,i,ilev) = rand_num_mt - enddo - enddo - enddo - endif - - do ilev = 2,nlay - do i = 1, ncol - do isubcol = 1, nsubcol - if (CDF(isubcol, i, ilev-1) > 1._rb - cldf(i,ilev-1) )& - & then - CDF(isubcol,i,ilev) = CDF(isubcol,i,ilev-1) - else - CDF(isubcol,i,ilev) = CDF(isubcol,i,ilev) * (1._rb & - & - cldf(i,ilev-1)) - endif - enddo - enddo - enddo - - case(3) -! Maximum overlap -! i) pick the same random numebr at every level - - if (irng.eq.0) then - do isubcol = 1,nsubcol - call kissvec(seed1, seed2, seed3, seed4, rand_num) - do ilev = 1,nlay - CDF(isubcol,:,ilev) = rand_num - enddo - enddo - elseif (irng.eq.1) then - do isubcol = 1, nsubcol - do i = 1, ncol - rand_num_mt = getRandomReal(randomNumbers) - do ilev = 1, nlay - CDF(isubcol,i,ilev) = rand_num_mt - enddo - enddo - enddo - endif - -! mji - Activate exponential cloud overlap option - case(4) - ! Exponential overlap: weighting between maximum and random overlap increases with the distance. - ! The random numbers for exponential overlap verify: - ! j=1 RAN(j)=RND1 - ! j>1 if RND1 < alpha(j,j-1) => RAN(j) = RAN(j-1) - ! RAN(j) = RND2 - ! alpha is obtained from the equation - ! alpha = exp(-(Z(j)-Z(j-1))/Zo) where Zo is a characteristic length scale - - ! compute alpha - do i = 1, ncol - alpha(i, 1) = 0._rb - do ilev = 2,nlay - alpha(i, ilev) = exp( -( hgt (i, ilev) - & - & hgt (i, ilev-1)) / Zo) - enddo - enddo - - ! generate 2 streams of random numbers - if (irng.eq.0) then - do isubcol = 1,nsubcol - do ilev = 1,nlay - call kissvec(seed1, seed2, seed3, seed4, rand_num) - CDF(isubcol, :, ilev) = rand_num - call kissvec(seed1, seed2, seed3, seed4, rand_num) - CDF2(isubcol, :, ilev) = rand_num - enddo - enddo - elseif (irng.eq.1) then - do isubcol = 1, nsubcol - do i = 1, ncol - do ilev = 1, nlay - rand_num_mt = getRandomReal(randomNumbers) - CDF(isubcol,i,ilev) = rand_num_mt - rand_num_mt = getRandomReal(randomNumbers) - CDF2(isubcol,i,ilev) = rand_num_mt - enddo - enddo - enddo - endif - - ! generate random numbers - do ilev = 2,nlay - where (CDF2(:, :, ilev) < spread(alpha (:,ilev), & - & dim=1,nCopies=nsubcol) ) - CDF(:,:,ilev) = CDF(:,:,ilev-1) - end where - end do - -! Activate exponential-random cloud overlap option - case(5) - ! Exponential-random overlap: -!mz* call wrf_error_fatal("Cloud Overlap case 5: ER has not yet & -! been implemented. Stopping...") - - end select - -! -- generate subcolumns for homogeneous clouds ----- - do ilev = 1,nlay - iscloudy(:,:,ilev) = (CDF(:,:,ilev) >= 1._rb - & - & spread(cldf(:,ilev), dim=1, nCopies=nsubcol) ) - enddo - -! where the subcolumn is cloudy, the subcolumn cloud fraction is 1; -! where the subcolumn is not cloudy, the subcolumn cloud fraction is 0; -! where there is a cloud, define the subcolumn cloud properties, -! otherwise set these to zero - - do ilev = 1,nlay - do i = 1, ncol - do isubcol = 1, nsubcol - if (iscloudy(isubcol,i,ilev) ) then - cld_stoch(isubcol,i,ilev) = 1._rb - clwp_stoch(isubcol,i,ilev) = clwp(i,ilev) - ciwp_stoch(isubcol,i,ilev) = ciwp(i,ilev) -!mz -! cswp_stoch(isubcol,i,ilev) = cswp(i,ilev) - cswp_stoch(isubcol,i,ilev) = 0._rb - n = ngb(isubcol) - tauc_stoch(isubcol,i,ilev) = tauc(n,i,ilev) -! ssac_stoch(isubcol,i,ilev) = ssac(n,i,ilev) -! asmc_stoch(isubcol,i,ilev) = asmc(n,i,ilev) - else - cld_stoch(isubcol,i,ilev) = 0._rb - clwp_stoch(isubcol,i,ilev) = 0._rb - ciwp_stoch(isubcol,i,ilev) = 0._rb - cswp_stoch(isubcol,i,ilev) = 0._rb - tauc_stoch(isubcol,i,ilev) = 0._rb -! ssac_stoch(isubcol,i,ilev) = 1._rb -! asmc_stoch(isubcol,i,ilev) = 1._rb - endif - enddo - enddo - enddo - -! -- compute the means of the subcolumns --- -! mean_cld_stoch(:,:) = 0._rb -! mean_clwp_stoch(:,:) = 0._rb -! mean_ciwp_stoch(:,:) = 0._rb -! mean_tauc_stoch(:,:) = 0._rb -! mean_ssac_stoch(:,:) = 0._rb -! mean_asmc_stoch(:,:) = 0._rb -! do i = 1, nsubcol -! mean_cld_stoch(:,:) = cld_stoch(i,:,:) + mean_cld_stoch(:,:) -! mean_clwp_stoch(:,:) = clwp_stoch( i,:,:) + mean_clwp_stoch(:,:) -! mean_ciwp_stoch(:,:) = ciwp_stoch( i,:,:) + mean_ciwp_stoch(:,:) -! mean_tauc_stoch(:,:) = tauc_stoch( i,:,:) + mean_tauc_stoch(:,:) -! mean_ssac_stoch(:,:) = ssac_stoch( i,:,:) + mean_ssac_stoch(:,:) -! mean_asmc_stoch(:,:) = asmc_stoch( i,:,:) + mean_asmc_stoch(:,:) -! end do -! mean_cld_stoch(:,:) = mean_cld_stoch(:,:) / nsubcol -! mean_clwp_stoch(:,:) = mean_clwp_stoch(:,:) / nsubcol -! mean_ciwp_stoch(:,:) = mean_ciwp_stoch(:,:) / nsubcol -! mean_tauc_stoch(:,:) = mean_tauc_stoch(:,:) / nsubcol -! mean_ssac_stoch(:,:) = mean_ssac_stoch(:,:) / nsubcol -! mean_asmc_stoch(:,:) = mean_asmc_stoch(:,:) / nsubcol - - end subroutine generate_stochastic_clouds - -!------------------------------------------------------------------ -! Private subroutines -!------------------------------------------------------------------ - -!----------------------------------------------------------------- - subroutine kissvec(seed1,seed2,seed3,seed4,ran_arr) -!---------------------------------------------------------------- - -! public domain code -! made available from http://www.fortran.com/ -! downloaded by pjr on 03/16/04 for NCAR CAM -! converted to vector form, functions inlined by pjr,mvr on 05/10/2004 - -! The KISS (Keep It Simple Stupid) random number generator. Combines: -! (1) The congruential generator x(n)=69069*x(n-1)+1327217885, period 2^32. -! (2) A 3-shift shift-register generator, period 2^32-1, -! (3) Two 16-bit multiply-with-carry generators, period 597273182964842497>2^59 -! Overall period>2^123; - real(kind=rb), dimension(:), intent(inout) :: ran_arr - integer(kind=im), dimension(:), intent(inout) :: seed1,seed2,seed3& - & ,seed4 - integer(kind=im) :: i,sz,kiss - integer(kind=im) :: m, k, n - -! inline function - m(k, n) = ieor (k, ishft (k, n) ) - - sz = size(ran_arr) - do i = 1, sz - seed1(i) = 69069_im * seed1(i) + 1327217885_im - seed2(i) = m (m (m (seed2(i), 13_im), - 17_im), 5_im) - seed3(i) = 18000_im * iand (seed3(i), 65535_im) + & - & ishft (seed3(i), - 16_im) - seed4(i) = 30903_im * iand (seed4(i), 65535_im) + & - & ishft (seed4(i), - 16_im) - kiss = seed1(i) + seed2(i) + ishft (seed3(i), 16_im) + seed4(i) - ran_arr(i) = kiss*2.328306e-10_rb + 0.5_rb - end do - - end subroutine kissvec -! - subroutine rtrnmc_mcica(nlayers, istart, iend, iout, pz, semiss, & - & ncbands, cldfmc, taucmc, planklay, planklev, &!plankbnd, & - & pwvcm, fracs, taut, & - & totuflux, totdflux, htr, & - & totuclfl, totdclfl, htrc ) -!--------------------------------------------------------------- -! -! Original version: E. J. Mlawer, et al. RRTM_V3.0 -! Revision for GCMs: Michael J. Iacono; October, 2002 -! Revision for F90: Michael J. Iacono; June, 2006 -! -! This program calculates the upward fluxes, downward fluxes, and -! heating rates for an arbitrary clear or cloudy atmosphere. The input -! to this program is the atmospheric profile, all Planck function -! information, and the cloud fraction by layer. A variable diffusivity -! angle (SECDIFF) is used for the angle integration. Bands 2-3 and 5-9 -! use a value for SECDIFF that varies from 1.50 to 1.80 as a function of -! the column water vapor, and other bands use a value of 1.66. The Gaussian -! weight appropriate to this angle (WTDIFF=0.5) is applied here. Note that -! use of the emissivity angle for the flux integration can cause errors of -! 1 to 4 W/m2 within cloudy layers. -! Clouds are treated with the McICA stochastic approach and maximum-random -! cloud overlap. -!*************************************************************************** - -! ------- Declarations ------- - -! ----- Input ----- - integer(kind=im), intent(in) :: nlayers ! total number of layers - integer(kind=im), intent(in) :: istart ! beginning band of calculation - integer(kind=im), intent(in) :: iend ! ending band of calculation - integer(kind=im), intent(in) :: iout ! output option flag - -! Atmosphere - real(kind=rb), intent(in) :: pz(0:) ! level (interface) pressures (hPa, mb) - ! Dimensions: (0:nlayers) - real(kind=rb), intent(in) :: pwvcm ! precipitable water vapor (cm) - real(kind=rb), intent(in) :: semiss(:) ! lw surface emissivity - ! Dimensions: (nbndlw) -!mz - real(kind=rb), intent(in) :: planklay(0:,:) ! - ! Dimensions: (nlayers,nbndlw) - real(kind=rb), intent(in) :: planklev(0:,:) ! - ! Dimensions: (0:nlayers,nbndlw) -! real(kind=rb), intent(in) :: plankbnd(:) ! - ! Dimensions: (nbndlw) - real(kind=rb), intent(in) :: fracs(:,:) ! - ! Dimensions: (nlayers,ngptw) - real(kind=rb), intent(in) :: taut(:,:) ! gaseous + aerosol optical depths - ! Dimensions: (nlayers,ngptlw) - -! Clouds - integer(kind=im), intent(in) :: ncbands ! number of cloud spectral bands - real(kind=rb), intent(in) :: cldfmc(:,:) ! layer cloud fraction [mcica] - ! Dimensions: (ngptlw,nlayers) - real(kind=rb), intent(in) :: taucmc(:,:) ! layer cloud optical depth [mcica] - ! Dimensions: (ngptlw,nlayers) - -! ----- Output ----- - real(kind=rb), intent(out) :: totuflux(0:) ! upward longwave flux (w/m2) - ! Dimensions: (0:nlayers) - real(kind=rb), intent(out) :: totdflux(0:) ! downward longwave flux (w/m2) - ! Dimensions: (0:nlayers) -!mz* real(kind=rb), intent(out) :: fnet(0:) ! net longwave flux (w/m2) - ! Dimensions: (0:nlayers) - real(kind=rb), intent(out) :: htr(:) -!mz real(kind=rb), intent(out) :: htr(0:) ! longwave heating rate (k/day) - ! Dimensions: (0:nlayers) - real(kind=rb), intent(out) :: totuclfl(0:) ! clear sky upward longwave flux (w/m2) - ! Dimensions: (0:nlayers) - real(kind=rb), intent(out) :: totdclfl(0:) ! clear sky downward longwave flux (w/m2) - ! Dimensions: (0:nlayers) -!mz*real(kind=rb), intent(out) :: fnetc(0:) ! clear sky net longwave flux (w/m2) - ! Dimensions: (0:nlayers) - real(kind=rb), intent(out) :: htrc(:) -! real(kind=rb), intent(out) :: htrc(0:) ! clear sky longwave heating rate (k/day) - ! Dimensions: (0:nlayers) - -! ----- Local ----- -! Declarations for radiative transfer - real (kind=kind_phys), dimension(0:nlayers) :: fnet, fnetc - real(kind=rb) :: abscld(nlayers,ngptlw) - real(kind=rb) :: atot(nlayers) - real(kind=rb) :: atrans(nlayers) - real(kind=rb) :: bbugas(nlayers) - real(kind=rb) :: bbutot(nlayers) - real(kind=rb) :: clrurad(0:nlayers) - real(kind=rb) :: clrdrad(0:nlayers) - real(kind=rb) :: efclfrac(nlayers,ngptlw) - real(kind=rb) :: uflux(0:nlayers) - real(kind=rb) :: dflux(0:nlayers) - real(kind=rb) :: urad(0:nlayers) - real(kind=rb) :: drad(0:nlayers) - real(kind=rb) :: uclfl(0:nlayers) - real(kind=rb) :: dclfl(0:nlayers) - real(kind=rb) :: odcld(nlayers,ngptlw) - - - real(kind=rb) :: secdiff(nbands) ! secant of diffusivity angle - real(kind=rb) :: transcld, radld, radclrd, plfrac, blay, dplankup,& - & dplankdn - real(kind=rb) :: odepth, odtot, odepth_rec, odtot_rec, gassrc - real(kind=rb) :: tblind, tfactot, bbd, bbdtot, tfacgas, transc, & - & tausfac - real(kind=rb) :: rad0, reflect, radlu, radclru - - integer(kind=im) :: icldlyr(nlayers) ! flag for cloud in layer - integer(kind=im) :: ibnd, ib, iband, lay, lev, l, ig ! loop indices - integer(kind=im) :: igc ! g-point interval counter - integer(kind=im) :: iclddn ! flag for cloud in down path - integer(kind=im) :: ittot, itgas, itr ! lookup table indices -!mz* - real (kind=kind_phys), parameter :: rec_6 = 0.166667 - ! The cumulative sum of new g-points for each band - integer(kind=im) :: ngs(nbands) - ngs(:) = (/10,22,38,52,68,76,88,96,108,114,122,130,134,136,138, & - & 140/) - -! ------- Definitions ------- -! input -! nlayers ! number of model layers -! ngptlw ! total number of g-point subintervals -! nbndlw ! number of longwave spectral bands -! ncbands ! number of spectral bands for clouds -! secdiff ! diffusivity angle -! wtdiff ! weight for radiance to flux conversion -! pavel ! layer pressures (mb) -! pz ! level (interface) pressures (mb) -! tavel ! layer temperatures (k) -! tz ! level (interface) temperatures(mb) -! tbound ! surface temperature (k) -! cldfrac ! layer cloud fraction -! taucloud ! layer cloud optical depth -! itr ! integer look-up table index -! icldlyr ! flag for cloudy layers -! iclddn ! flag for cloud in column at any layer -! semiss ! surface emissivities for each band -! reflect ! surface reflectance -! bpade ! 1/(pade constant) -! tau_tbl ! clear sky optical depth look-up table -! exp_tbl ! exponential look-up table for transmittance -! tfn_tbl ! tau transition function look-up table - -! local -! atrans ! gaseous absorptivity -! abscld ! cloud absorptivity -! atot ! combined gaseous and cloud absorptivity -! odclr ! clear sky (gaseous) optical depth -! odcld ! cloud optical depth -! odtot ! optical depth of gas and cloud -! tfacgas ! gas-only pade factor, used for planck fn -! tfactot ! gas and cloud pade factor, used for planck fn -! bbdgas ! gas-only planck function for downward rt -! bbugas ! gas-only planck function for upward rt -! bbdtot ! gas and cloud planck function for downward rt -! bbutot ! gas and cloud planck function for upward calc. -! gassrc ! source radiance due to gas only -! efclfrac ! effective cloud fraction -! radlu ! spectrally summed upward radiance -! radclru ! spectrally summed clear sky upward radiance -! urad ! upward radiance by layer -! clrurad ! clear sky upward radiance by layer -! radld ! spectrally summed downward radiance -! radclrd ! spectrally summed clear sky downward radiance -! drad ! downward radiance by layer -! clrdrad ! clear sky downward radiance by layer - - -! output -! totuflux ! upward longwave flux (w/m2) -! totdflux ! downward longwave flux (w/m2) -! fnet ! net longwave flux (w/m2) -! htr ! longwave heating rate (k/day) -! totuclfl ! clear sky upward longwave flux (w/m2) -! totdclfl ! clear sky downward longwave flux (w/m2) -! fnetc ! clear sky net longwave flux (w/m2) -! htrc ! clear sky longwave heating rate (k/day) - - -!jm not thread safe hvrrtc = '$Revision: 1.3 $' - - do ibnd = 1,nbands!mz*nbndlw - if (ibnd.eq.1 .or. ibnd.eq.4 .or. ibnd.ge.10) then - secdiff(ibnd) = 1.66_rb - else - secdiff(ibnd) = a0(ibnd) + a1(ibnd)*exp(a2(ibnd)*pwvcm) - if (secdiff(ibnd) .gt. 1.80_rb) secdiff(ibnd) = 1.80_rb - if (secdiff(ibnd) .lt. 1.50_rb) secdiff(ibnd) = 1.50_rb - endif - enddo - - urad(0) = 0.0_rb - drad(0) = 0.0_rb - totuflux(0) = 0.0_rb - totdflux(0) = 0.0_rb - clrurad(0) = 0.0_rb - clrdrad(0) = 0.0_rb - totuclfl(0) = 0.0_rb - totdclfl(0) = 0.0_rb - - do lay = 1, nlayers - urad(lay) = 0.0_rb - drad(lay) = 0.0_rb - totuflux(lay) = 0.0_rb - totdflux(lay) = 0.0_rb - clrurad(lay) = 0.0_rb - clrdrad(lay) = 0.0_rb - totuclfl(lay) = 0.0_rb - totdclfl(lay) = 0.0_rb - icldlyr(lay) = 0 - -! Change to band loop? - do ig = 1, ngptlw - if (cldfmc(ig,lay) .eq. 1._rb) then - ib = ngb(ig) - odcld(lay,ig) = secdiff(ib) * taucmc(ig,lay) - transcld = exp(-odcld(lay,ig)) - abscld(lay,ig) = 1._rb - transcld - efclfrac(lay,ig) = abscld(lay,ig) * cldfmc(ig,lay) - icldlyr(lay) = 1 - else - odcld(lay,ig) = 0.0_rb - abscld(lay,ig) = 0.0_rb - efclfrac(lay,ig) = 0.0_rb - endif - enddo - - enddo - - igc = 1 -! Loop over frequency bands. - do iband = istart, iend - -! Reinitialize g-point counter for each band if output for each band is requested. - if (iout.gt.0.and.iband.ge.2) igc = ngs(iband-1)+1 - -! Loop over g-channels. - 1000 continue - -! Radiative transfer starts here. - radld = 0._rb - radclrd = 0._rb - iclddn = 0 - -! Downward radiative transfer loop. - - do lev = nlayers, 1, -1 - plfrac = fracs(lev,igc) - blay = planklay(lev,iband) - dplankup = planklev(lev,iband) - blay - dplankdn = planklev(lev-1,iband) - blay - odepth = secdiff(iband) * taut(lev,igc) - if (odepth .lt. 0.0_rb) odepth = 0.0_rb -! Cloudy layer - if (icldlyr(lev).eq.1) then - iclddn = 1 - odtot = odepth + odcld(lev,igc) - if (odtot .lt. 0.06_rb) then - atrans(lev) = odepth - 0.5_rb*odepth*odepth - odepth_rec = rec_6*odepth - gassrc = plfrac*(blay+dplankdn*odepth_rec)*atrans(lev) - - atot(lev) = odtot - 0.5_rb*odtot*odtot - odtot_rec = rec_6*odtot - bbdtot = plfrac * (blay+dplankdn*odtot_rec) - bbd = plfrac*(blay+dplankdn*odepth_rec) - radld = radld - radld * (atrans(lev) + & - & efclfrac(lev,igc) * (1. - atrans(lev))) + & - & gassrc + cldfmc(igc,lev) * & - & (bbdtot * atot(lev) - gassrc) - drad(lev-1) = drad(lev-1) + radld - - bbugas(lev) = plfrac * (blay+dplankup*odepth_rec) - bbutot(lev) = plfrac * (blay+dplankup*odtot_rec) - - elseif (odepth .le. 0.06_rb) then - atrans(lev) = odepth - 0.5_rb*odepth*odepth - odepth_rec = rec_6*odepth - gassrc = plfrac*(blay+dplankdn*odepth_rec)*atrans(lev) - - odtot = odepth + odcld(lev,igc) - tblind = odtot/(bpade+odtot) - ittot = tblint*tblind + 0.5_rb - tfactot = tfn_tbl(ittot) - bbdtot = plfrac * (blay + tfactot*dplankdn) - bbd = plfrac*(blay+dplankdn*odepth_rec) - atot(lev) = 1. - exp_tbl(ittot) - - radld = radld - radld * (atrans(lev) + & - & efclfrac(lev,igc) * (1._rb - atrans(lev))) + & - & gassrc + cldfmc(igc,lev) * & - & (bbdtot * atot(lev) - gassrc) - drad(lev-1) = drad(lev-1) + radld - - bbugas(lev) = plfrac * (blay + dplankup*odepth_rec) - bbutot(lev) = plfrac * (blay + tfactot * dplankup) - - else - - tblind = odepth/(bpade+odepth) - itgas = tblint*tblind+0.5_rb - odepth = tau_tbl(itgas) - atrans(lev) = 1._rb - exp_tbl(itgas) - tfacgas = tfn_tbl(itgas) - gassrc = atrans(lev) * plfrac * (blay + tfacgas*dplankdn) - - odtot = odepth + odcld(lev,igc) - tblind = odtot/(bpade+odtot) - ittot = tblint*tblind + 0.5_rb - tfactot = tfn_tbl(ittot) - bbdtot = plfrac * (blay + tfactot*dplankdn) - bbd = plfrac*(blay+tfacgas*dplankdn) - atot(lev) = 1._rb - exp_tbl(ittot) - - radld = radld - radld * (atrans(lev) + & - & efclfrac(lev,igc) * (1._rb - atrans(lev))) + & - & gassrc + cldfmc(igc,lev) * & - & (bbdtot * atot(lev) - gassrc) - drad(lev-1) = drad(lev-1) + radld - bbugas(lev) = plfrac * (blay + tfacgas * dplankup) - bbutot(lev) = plfrac * (blay + tfactot * dplankup) - endif -! Clear layer - else - if (odepth .le. 0.06_rb) then - atrans(lev) = odepth-0.5_rb*odepth*odepth - odepth = rec_6*odepth - bbd = plfrac*(blay+dplankdn*odepth) - bbugas(lev) = plfrac*(blay+dplankup*odepth) - else - tblind = odepth/(bpade+odepth) - itr = tblint*tblind+0.5_rb - transc = exp_tbl(itr) - atrans(lev) = 1._rb-transc - tausfac = tfn_tbl(itr) - bbd = plfrac*(blay+tausfac*dplankdn) - bbugas(lev) = plfrac * (blay + tausfac * dplankup) - endif - radld = radld + (bbd-radld)*atrans(lev) - drad(lev-1) = drad(lev-1) + radld - endif -! Set clear sky stream to total sky stream as long as layers -! remain clear. Streams diverge when a cloud is reached (iclddn=1), -! and clear sky stream must be computed separately from that point. - if (iclddn.eq.1) then - radclrd = radclrd + (bbd-radclrd) * atrans(lev) - clrdrad(lev-1) = clrdrad(lev-1) + radclrd - else - radclrd = radld - clrdrad(lev-1) = drad(lev-1) - endif - enddo - -! Spectral emissivity & reflectance -! Include the contribution of spectrally varying longwave emissivity -! and reflection from the surface to the upward radiative transfer. -! Note: Spectral and Lambertian reflection are identical for the -! diffusivity angle flux integration used here. - -!mz* -! rad0 = fracs(1,igc) * plankbnd(iband) - rad0 = semiss(iband) * fracs(1,igc) * planklay(0,iband) -!mz -! Add in specular reflection of surface downward radiance. - reflect = 1._rb - semiss(iband) - radlu = rad0 + reflect * radld - radclru = rad0 + reflect * radclrd - - -! Upward radiative transfer loop. - urad(0) = urad(0) + radlu - clrurad(0) = clrurad(0) + radclru - - do lev = 1, nlayers -! Cloudy layer - if (icldlyr(lev) .eq. 1) then - gassrc = bbugas(lev) * atrans(lev) - radlu = radlu - radlu * (atrans(lev) + & - & efclfrac(lev,igc) * (1._rb - atrans(lev))) + & - & gassrc + cldfmc(igc,lev) * & - & (bbutot(lev) * atot(lev) - gassrc) - urad(lev) = urad(lev) + radlu -! Clear layer - else - radlu = radlu + (bbugas(lev)-radlu)*atrans(lev) - urad(lev) = urad(lev) + radlu - endif -! Set clear sky stream to total sky stream as long as all layers -! are clear (iclddn=0). Streams must be calculated separately at -! all layers when a cloud is present (ICLDDN=1), because surface -! reflectance is different for each stream. - if (iclddn.eq.1) then - radclru = radclru + (bbugas(lev)-radclru)*atrans(lev) - clrurad(lev) = clrurad(lev) + radclru - else - radclru = radlu - clrurad(lev) = urad(lev) - endif - enddo - -! Increment g-point counter - igc = igc + 1 -! Return to continue radiative transfer for all g-channels in present band - if (igc .le. ngs(iband)) go to 1000 - -! Process longwave output from band for total and clear streams. -! Calculate upward, downward, and net flux. - do lev = nlayers, 0, -1 - uflux(lev) = urad(lev)*wtdiff - dflux(lev) = drad(lev)*wtdiff - urad(lev) = 0.0_rb - drad(lev) = 0.0_rb - totuflux(lev) = totuflux(lev) + uflux(lev) * delwave(iband) - totdflux(lev) = totdflux(lev) + dflux(lev) * delwave(iband) - uclfl(lev) = clrurad(lev)*wtdiff - dclfl(lev) = clrdrad(lev)*wtdiff - clrurad(lev) = 0.0_rb - clrdrad(lev) = 0.0_rb - totuclfl(lev) = totuclfl(lev) + uclfl(lev) * delwave(iband) - totdclfl(lev) = totdclfl(lev) + dclfl(lev) * delwave(iband) - enddo - -! End spectral band loop - enddo - -! Calculate fluxes at surface - totuflux(0) = totuflux(0) * fluxfac - totdflux(0) = totdflux(0) * fluxfac - fnet(0) = totuflux(0) - totdflux(0) - totuclfl(0) = totuclfl(0) * fluxfac - totdclfl(0) = totdclfl(0) * fluxfac - fnetc(0) = totuclfl(0) - totdclfl(0) - -! Calculate fluxes at model levels - do lev = 1, nlayers - totuflux(lev) = totuflux(lev) * fluxfac - totdflux(lev) = totdflux(lev) * fluxfac - fnet(lev) = totuflux(lev) - totdflux(lev) - totuclfl(lev) = totuclfl(lev) * fluxfac - totdclfl(lev) = totdclfl(lev) * fluxfac - fnetc(lev) = totuclfl(lev) - totdclfl(lev) - l = lev - 1 - -! Calculate heating rates at model layers - htr(l)=heatfac*(fnet(l)-fnet(lev))/(pz(l)-pz(lev)) - htrc(l)=heatfac*(fnetc(l)-fnetc(lev))/(pz(l)-pz(lev)) - enddo - -! Set heating rate to zero in top layer - htr(nlayers) = 0.0_rb - htrc(nlayers) = 0.0_rb - - end subroutine rtrnmc_mcica ! ------------------------------------------------------------------------------ subroutine cldprmc(nlayers, inflag, iceflag, liqflag, cldfmc, & @@ -8968,4 +7792,4 @@ end subroutine cldprmc !........................................!$ end module rrtmg_lw !$ -!========================================!$ \ No newline at end of file +!========================================!$ diff --git a/physics/radlw_main.meta b/physics/radlw_main.meta index 70d820922..d96ab9af1 100644 --- a/physics/radlw_main.meta +++ b/physics/radlw_main.meta @@ -1,7 +1,7 @@ [ccpp-table-properties] name = rrtmg_lw type = scheme - dependencies = mersenne_twister.f,physcons.F90,physparam.f,radlw_datatb.f,radlw_param.f,HWRF_mcica_random_numbers.F90,HWRF_mersenne_twister.F90 + dependencies = mersenne_twister.f,physcons.F90,physparam.f,radlw_datatb.f,radlw_param.f ######################################################################## [ccpp-arg-table] diff --git a/physics/radsw_main.F90 b/physics/radsw_main.F90 index 0f5a8b110..44de9848c 100644 --- a/physics/radsw_main.F90 +++ b/physics/radsw_main.F90 @@ -404,9 +404,7 @@ module rrtmg_sw ! --- public accessable subprograms - public rrtmg_sw_init, rrtmg_sw_run, rrtmg_sw_finalize, rswinit, & - & kissvec, generate_stochastic_clouds_sw, mcica_subcol_sw - + public rrtmg_sw_init, rrtmg_sw_run, rrtmg_sw_finalize, rswinit ! ================= contains @@ -758,74 +756,6 @@ subroutine rrtmg_sw_run & & intent(inout) :: fdncmp ! --- locals: -!mz* HWRF -- input of mcica_subcol_sw - real(kind=kind_phys),dimension(npts,nlay) :: hgt - real(kind=kind_phys) :: dzsum - real(kind=kind_phys),dimension( nbdsw, npts, nlay ) :: taucld3, & - ssacld3, & - asmcld3, & - fsfcld3 - -!mz* HWRF -- OUTPUT from mcica_subcol_sw - real(kind=kind_phys),dimension(ngptsw,npts,nlay) :: cldfmcl ! Cloud fraction - ! Dimensions: (ngptsw,ncol,nlay) - real(kind=kind_phys),dimension(ngptsw,npts,nlay) :: ciwpmcl ! In-cloud ice water path (g/m2) - ! Dimensions: (ngptsw,ncol,nlay) - real(kind=kind_phys),dimension(ngptsw,npts,nlay) :: clwpmcl ! In-cloud liquid water path (g/m2) - ! Dimensions: (ngptsw,ncol,nlay) - real(kind=kind_phys),dimension(ngptsw,npts,nlay) :: cswpmcl ! In-cloud snow water path (g/m2) - ! Dimensions: (ngptsw,ncol,nlay) - real(kind=kind_phys),dimension(npts,nlay) :: relqmcl ! Cloud water drop effective radius (microns) - ! Dimensions: (ncol,nlay) - real(kind=kind_phys),dimension(npts,nlay) :: reicmcl ! Cloud ice effective size (microns) - ! Dimensions: (ncol,nlay) - real(kind=kind_phys),dimension(npts,nlay) :: resnmcl ! Snow effective size (microns) - ! Dimensions: (ncol,nlay) - real(kind=kind_phys),dimension(ngptsw,npts,nlay) :: taucmcl ! In-cloud optical depth - ! Dimensions: (ngptsw,ncol,nlay) - real(kind=kind_phys),dimension(ngptsw,npts,nlay) :: ssacmcl ! in-cloud single scattering albedo [mcica] - ! Dimensions: (ngptsw,ncol,nlay) - real(kind=kind_phys),dimension(ngptsw,npts,nlay) :: asmcmcl ! in-cloud asymmetry parameter [mcica] - ! Dimensions: (ngptsw,ncol,nlay) - real(kind=kind_phys),dimension(ngptsw,npts,nlay) :: fsfcmcl ! in-cloud forward scattering fraction [mcica] - ! Dimensions: (ngptsw,ncol,nlay) -!HWRF cldprmc_sw input -! real(kind=kind_phys),dimension(ngptsw,nlay) :: cldfmc,cldfmc_save! cloud fraction [mcica] -! ! Dimensions: (ngptsw,nlayers) - real(kind=kind_phys),dimension(ngptsw,nlay) :: ciwpmc ! cloud ice water path [mcica] - ! Dimensions: (ngptsw,nlayers) - real(kind=kind_phys),dimension(ngptsw,nlay) :: clwpmc ! cloud liquid water path [mcica] - ! Dimensions: (ngptsw,nlayers) - real(kind=kind_phys),dimension(ngptsw,nlay) :: cswpmc ! cloud snow water path [mcica] - ! Dimensions: (ngptsw,nlayers) - real(kind=kind_phys),dimension(nlay) :: resnmc ! cloud snow particle effective radius (microns) - ! Dimensions: (nlayers) - real(kind=kind_phys),dimension(nlay) :: relqmc ! cloud liquid particle effective radius (microns) - ! Dimensions: (nlayers) - real(kind=kind_phys),dimension(nlay) :: reicmc ! cloud ice particle effective radius (microns) - ! Dimensions: (nlayers) - ! specific definition of reicmc depends on setting of iceflag: - ! iceflag = 1: ice effective radius, r_ec, (Ebert and Curry, 1992), - ! r_ec range is limited to 13.0 to 130.0 microns - ! iceflag = 2: ice effective radius, r_k, (Key, Streamer Ref. Manual, 1996) - ! r_k range is limited to 5.0 to 131.0 microns - ! iceflag = 3: generalized effective size, dge, (Fu, 1996), - ! dge range is limited to 5.0 to 140.0 microns - ! [dge = 1.0315 * r_ec] - real(kind=kind_phys),dimension(ngptsw,nlay) :: fsfcmc ! cloud forward scattering fraction - ! Dimensions: (ngptsw,nlayers) - -!mz* HWRF cldprmc_sw output (delta scaled) - real(kind=kind_phys),dimension(ngptsw,nlay) :: taucmc ! cloud optical depth (delta scaled) - ! Dimensions: (ngptsw,nlayers) - real(kind=kind_phys),dimension(ngptsw,nlay) :: ssacmc ! single scattering albedo (delta scaled) - ! Dimensions: (ngptsw,nlayers) - real(kind=kind_phys),dimension(ngptsw,nlay) :: asmcmc ! asymmetry parameter (delta scaled) - ! Dimensions: (ngptsw,nlayers) - real(kind=kind_phys),dimension(ngptsw,nlay) :: taormc ! cloud optical depth (non-delta scaled) - ! Dimensions: (ngptsw,nlayers) -!mz* - real (kind=kind_phys), dimension(nlay,ngptsw) :: cldfmc, & & cldfmc_save, & & taug, taur @@ -980,52 +910,6 @@ subroutine rrtmg_sw_run & albbm(2) = sfcalb_uvis_dir(j1) albdf(2) = sfcalb_uvis_dif(j1) - -! mz*: HWRF - if (iovr == 4 ) then - - -!Add layer height needed for exponential (icld=4) and -! exponential-random (icld=5) overlap options - - !iplon = 1 - irng = 0 - permuteseed = 1 - -!mz* Derive height of each layer mid-point from layer thickness. -! Needed for exponential (iovr=4) and exponential-random overlap -! option (iovr=5)only. - dzsum =0.0 - do k = 1,nlay - hgt(j1,k)= dzsum+0.5*dzlyr(j1,k)*1000. !km->m - dzsum = dzsum+ dzlyr(j1,k)*1000. - enddo - -! Zero out cloud optical properties here; not used when passing physical properties -! to radiation and taucld is calculated in radiation - do k = 1, nlay - do ib = 1, nbdsw - taucld3(ib,j1,k) = 0.0 - ssacld3(ib,j1,k) = 1.0 - asmcld3(ib,j1,k) = 0.0 - fsfcld3(ib,j1,k) = 0.0 - enddo - enddo - - call mcica_subcol_sw (1, 1, nlay, iovr, permuteseed, & - & irng, plyr(j1:j1,:), hgt(j1:j1,:), & - & cld_cf(j1:j1,:), cld_iwp(j1:j1,:), cld_lwp(j1:j1,:), & - & cld_swp(j1:j1,:), cld_ref_ice(j1:j1,:), cld_ref_liq(j1:j1,:), & - & cld_ref_snow(j1:j1,:), taucld3(:,j1:j1,:), ssacld3(:,j1:j1,:), & - & asmcld3(:,j1:j1,:), fsfcld3(:,j1:j1,:), cldfmcl(:,j1:j1,:), & !--output - & ciwpmcl(:,j1:j1,:), clwpmcl(:,j1:j1,:), cswpmcl(:,j1:j1,:), & - & reicmcl(j1:j1,:), relqmcl(j1:j1,:), resnmcl(j1:j1,:), & - & taucmcl(:,j1:j1,:), ssacmcl(:,j1:j1,:), asmcmcl(:,j1:j1,:), & - & fsfcmcl(:,j1:j1,:)) - - endif -!mz* end - !> -# Prepare atmospheric profile for use in rrtm. ! the vertical index of internal array is from surface to top @@ -1111,28 +995,6 @@ subroutine rrtmg_sw_run & cdat3(k) = cld_swp(j1,kk) ! cloud snow path cdat4(k) = cld_ref_snow(j1,kk) ! snow partical effctive radius enddo - if (iovr == 4) then !mz* HWRF - do k = 1, nlay - kk = nlp1 - k - do ig = 1, ngptsw - cldfmc(k,ig) = cldfmcl(ig,j1,kk) - taucmc(ig,k) = taucmcl(ig,j1,kk) - ssacmc(ig,k) = ssacmcl(ig,j1,kk) - asmcmc(ig,k) = asmcmcl(ig,j1,kk) - fsfcmc(ig,k) = fsfcmcl(ig,j1,kk) - ciwpmc(ig,k) = ciwpmcl(ig,j1,kk) - clwpmc(ig,k) = clwpmcl(ig,j1,kk) - if (iceflgsw.eq.5) then - cswpmc(ig,k) = cswpmcl(ig,j1,kk) - endif - enddo - reicmc(k) = reicmcl(j1,kk) - relqmc(k) = relqmcl(j1,kk) - if (iceflgsw.eq.5) then - resnmc(k) = resnmcl(j1,kk) - endif - enddo - endif else ! use diagnostic cloud method do k = 1, nlay kk = nlp1 - k @@ -1226,31 +1088,6 @@ subroutine rrtmg_sw_run & cdat3(k) = cld_swp(j1,k) ! cloud snow path cdat4(k) = cld_ref_snow(j1,k) ! snow partical effctive radius enddo - if (iovr == 4) then !mz* HWRF -!mz* Move incoming GCM cloud arrays to RRTMG cloud arrays. -!For GCM input, incoming reicmcl is defined based on selected -!ice parameterization (inflglw) - do k = 1, nlay - do ig = 1, ngptsw - cldfmc(k,ig) = cldfmcl(ig,j1,k) - taucmc(ig,k) = taucmcl(ig,j1,k) - ssacmc(ig,k) = ssacmcl(ig,j1,k) - asmcmc(ig,k) = asmcmcl(ig,j1,k) - fsfcmc(ig,k) = fsfcmcl(ig,j1,k) - ciwpmc(ig,k) = ciwpmcl(ig,j1,k) - clwpmc(ig,k) = clwpmcl(ig,j1,k) - if (iceflgsw .eq. 5) then - cswpmc(ig,k) = cswpmcl(ig,j1,k) - endif - enddo - reicmc(k) = reicmcl(j1,k) - relqmc(k) = relqmcl(j1,k) - if (iceflgsw .eq. 5) then - resnmc(k) = resnmcl(j1,k) - endif - enddo - - end if else ! use diagnostic cloud method do k = 1, nlay cfrac(k) = cld_cf(j1,k) ! cloud fraction @@ -1273,7 +1110,7 @@ subroutine rrtmg_sw_run & do k = 1, nlay zcf0 = zcf0 * (f_one - cfrac(k)) enddo - else if (iovr == 1 .or. iovr == 4) then ! max/ran/exp overlapping + else if (iovr == 1) then ! max/ran/exp overlapping do k = 1, nlay if (cfrac(k) > ftiny) then ! cloudy layer zcf1 = min ( zcf1, f_one-cfrac(k) ) @@ -1298,15 +1135,6 @@ subroutine rrtmg_sw_run & if (zcf1 > f_zero) then ! cloudy sky column - !mz* for HWRF, save cldfmc with mcica - if (iovr == 4) then - do k = 1, nlay - do ig = 1, ngptsw - cldfmc_save(k,ig)=cldfmc (k,ig) - enddo - enddo - endif - call cldprop & ! --- inputs: & ( cfrac,cliqp,reliq,cicep,reice,cdat1,cdat2,cdat3,cdat4, & @@ -1315,15 +1143,6 @@ subroutine rrtmg_sw_run & & taucw, ssacw, asycw, cldfrc, cldfmc & & ) - if (iovr == 4) then - !mz for HWRF, still using mcica cldfmc - do k = 1, nlay - do ig = 1, ngptsw - cldfmc(k,ig)=cldfmc_save(k,ig) - enddo - enddo - endif - ! --- ... save computed layer cloud optical depth for output ! rrtm band 10 is approx to the 0.55 mu spectrum @@ -5804,578 +5623,6 @@ end subroutine taumol29 end subroutine taumol !----------------------------------- -!mz* HWRF subroutines - subroutine mcica_subcol_sw(iplon, ncol, nlay, icld, permuteseed, & - & irng, play, hgt, & - & cldfrac, ciwp, clwp, cswp, rei, rel, res, tauc, & - & ssac, asmc, fsfc, & - & cldfmcl, ciwpmcl, clwpmcl, cswpmcl, reicmcl, & - & relqmcl, resnmcl, & - & taucmcl, ssacmcl, asmcmcl, fsfcmcl) - -! ----- Input ----- -! Control - integer(kind=im), intent(in) :: iplon ! column/longitude dimension - integer(kind=im), intent(in) :: ncol ! number of columns - integer(kind=im), intent(in) :: nlay ! number of model layers - integer(kind=im), intent(in) :: icld ! clear/cloud, cloud overlap flag - integer(kind=im), intent(in) :: permuteseed ! if the cloud generator is called multiple times, - ! permute the seed between each call; - ! between calls for LW and SW, recommended - ! permuteseed differs by 'ngpt' - integer(kind=im), intent(inout) :: irng ! flag for random number generator - ! 0 = kissvec - ! 1 = Mersenne Twister - -! Atmosphere - real(kind=rb), intent(in) :: play(:,:) ! layer pressures (mb) - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: hgt(:,:) ! layer height (m) - ! Dimensions: (ncol,nlay) - -! Atmosphere/clouds - cldprop - real(kind=rb), intent(in) :: cldfrac(:,:) ! layer cloud fraction - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: tauc(:,:,:) ! in-cloud optical depth - ! Dimensions: (nbndsw,ncol,nlay) - real(kind=rb), intent(in) :: ssac(:,:,:) ! in-cloud single scattering albedo (non-delta scaled) - ! Dimensions: (nbndsw,ncol,nlay) - real(kind=rb), intent(in) :: asmc(:,:,:) ! in-cloud asymmetry parameter (non-delta scaled) - ! Dimensions: (nbndsw,ncol,nlay) - real(kind=rb), intent(in) :: fsfc(:,:,:) ! in-cloud forward scattering fraction (non-delta scaled) - ! Dimensions: (nbndsw,ncol,nlay) - real(kind=rb), intent(in) :: ciwp(:,:) ! in-cloud ice water path - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: clwp(:,:) ! in-cloud liquid water path - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: cswp(:,:) ! in-cloud snow water path - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: rei(:,:) ! cloud ice particle size - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: rel(:,:) ! cloud liquid particle size - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: res(:,:) ! cloud snow particle size - ! Dimensions: (ncol,nlay) - -! ----- Output ----- -! Atmosphere/clouds - cldprmc [mcica] - real(kind=rb), intent(out) :: cldfmcl(:,:,:) ! cloud fraction [mcica] - ! Dimensions: (ngptsw,ncol,nlay) - real(kind=rb), intent(out) :: ciwpmcl(:,:,:) ! in-cloud ice water path [mcica] - ! Dimensions: (ngptsw,ncol,nlay) - real(kind=rb), intent(out) :: clwpmcl(:,:,:) ! in-cloud liquid water path [mcica] - ! Dimensions: (ngptsw,ncol,nlay) - real(kind=rb), intent(out) :: cswpmcl(:,:,:) ! in-cloud snow water path [mcica] - ! Dimensions: (ngptsw,ncol,nlay) - real(kind=rb), intent(out) :: relqmcl(:,:) ! liquid particle size (microns) - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(out) :: reicmcl(:,:) ! ice partcle size (microns) - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(out) :: resnmcl(:,:) ! snow partcle size (microns) - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(out) :: taucmcl(:,:,:) ! in-cloud optical depth [mcica] - ! Dimensions: (ngptsw,ncol,nlay) - real(kind=rb), intent(out) :: ssacmcl(:,:,:) ! in-cloud single scattering albedo [mcica] - ! Dimensions: (ngptsw,ncol,nlay) - real(kind=rb), intent(out) :: asmcmcl(:,:,:) ! in-cloud asymmetry parameter [mcica] - ! Dimensions: (ngptsw,ncol,nlay) - real(kind=rb), intent(out) :: fsfcmcl(:,:,:) ! in-cloud forward scattering fraction [mcica] - ! Dimensions: (ngptsw,ncol,nlay) - -! ----- Local ----- - -! Stochastic cloud generator variables [mcica] - integer(kind=im), parameter :: nsubcsw = ngptsw ! number of sub-columns (g-point intervals) - integer(kind=im) :: ilev ! loop index - - real(kind=rb) :: pmid(ncol,nlay) ! layer pressures (Pa) -! real(kind=rb) :: pdel(ncol,nlay) ! layer pressure thickness (Pa) -! real(kind=rb) :: qi(ncol,nlay) ! ice water (specific humidity) -! real(kind=rb) :: ql(ncol,nlay) ! liq water (specific humidity) - -! Return if clear sky - if (icld.eq.0) return - -! NOTE: For GCM mode, permuteseed must be offset between LW and SW by at least number of subcolumns - -! Pass particle sizes to new arrays, no subcolumns for these properties yet -! Convert pressures from mb to Pa - - reicmcl(:ncol,:nlay) = rei(:ncol,:nlay) - relqmcl(:ncol,:nlay) = rel(:ncol,:nlay) - resnmcl(:ncol,:nlay) = res(:ncol,:nlay) - pmid(:ncol,:nlay) = play(:ncol,:nlay)*1.e2_rb - -! Convert input ice and liquid cloud water paths to specific humidity ice and liquid components - -! cwp = (q * pdel * 1000.) / gravit) -! = (kg/kg * kg m-1 s-2 *1000.) / m s-2 -! = (g m-2) -! -! q = (cwp * gravit) / (pdel *1000.) -! = (g m-2 * m s-2) / (kg m-1 s-2 * 1000.) -! = kg/kg - -! do ilev = 1, nlay -! qi(ilev) = (ciwp(ilev) * grav) / (pdel(ilev) * 1000._rb) -! ql(ilev) = (clwp(ilev) * grav) / (pdel(ilev) * 1000._rb) -! enddo - - call generate_stochastic_clouds_sw (ncol, nlay, nsubcsw, icld, & - & irng, pmid, hgt, cldfrac, clwp, ciwp, cswp, & - & tauc, ssac, asmc, fsfc, cldfmcl, clwpmcl, & - & ciwpmcl, cswpmcl, & - & taucmcl, ssacmcl, asmcmcl, fsfcmcl, permuteseed) - - end subroutine mcica_subcol_sw - -!------------------------------------------------------------------------------------------------- - subroutine generate_stochastic_clouds_sw(ncol, nlay, nsubcol, & - & icld, irng, pmid, hgt, cld, clwp, ciwp, cswp, & - & tauc, ssac, asmc, fsfc, cld_stoch, clwp_stoch, & - & ciwp_stoch, cswp_stoch, & - & tauc_stoch, ssac_stoch, asmc_stoch, fsfc_stoch, changeSeed) -!------------------------------------------------------------------------------------------------- -! Contact: Cecile Hannay (hannay@ucar.edu) -! -! Original code: Based on Raisanen et al., QJRMS, 2004. -! -! Modifications: Generalized for use with RRTMG and added Mersenne Twister as the default -! random number generator, which can be changed to the optional kissvec random number generator -! with flag 'irng'. Some extra functionality has been commented or removed. -! Michael J. Iacono, AER, Inc., February 2007 -! -! Given a profile of cloud fraction, cloud water and cloud ice, we produce a set of subcolumns. -! Each layer within each subcolumn is homogeneous, with cloud fraction equal to zero or one -! and uniform cloud liquid and cloud ice concentration. -! The ensemble as a whole reproduces the probability function of cloud liquid and ice within each layer -! and obeys an overlap assumption in the vertical. -! -! Overlap assumption: -! The cloud are consistent with 4 overlap assumptions: random, maximum, maximum-random and exponential. -! The default option is maximum-random (option 3) -! The options are: 1=random overlap, 2=max/random, 3=maximum overlap, 4=exponential overlap -! This is set with the variable "overlap" -!mji - Exponential overlap option (overlap=4) has been deactivated in this version -! The exponential overlap uses also a length scale, Zo. (real, parameter :: Zo = 2500. ) -! -! Seed: -! If the stochastic cloud generator is called several times during the same timestep, -! one should change the seed between the call to insure that the subcolumns are different. -! This is done by changing the argument 'changeSeed' -! For example, if one wants to create a set of columns for the shortwave and another set for the longwave , -! use 'changeSeed = 1' for the first call and'changeSeed = 2' for the second call -! -! PDF assumption: -! We can use arbitrary complicated PDFS. -! In the present version, we produce homogeneuous clouds (the simplest case). -! Future developments include using the PDF scheme of Ben Johnson. -! -! History file: -! Option to add diagnostics variables in the history file. (using FINCL in the namelist) -! nsubcol = number of subcolumns -! overlap = overlap type (1-3) -! Zo = length scale -! CLOUD_S = mean of the subcolumn cloud fraction ('_S" means Stochastic) -! CLDLIQ_S = mean of the subcolumn cloud water -! CLDICE_S = mean of the subcolumn cloud ice -! -! -! Note: -! Here: we force that the cloud condensate to be consistent with the cloud fraction -! i.e we only have cloud condensate when the cell is cloudy. -! In CAM: The cloud condensate and the cloud fraction are obtained from 2 different equations -! and the 2 quantities can be inconsistent (i.e. CAM can produce cloud fraction -! without cloud condensate or the opposite). -!---------------------------------------------------------------------- - - use mcica_random_numbers -! The Mersenne Twister random number engine - use MersenneTwister, only: randomNumberSequence, & - new_RandomNumberSequence, getRandomReal - - type(randomNumberSequence) :: randomNumbers - -! -- Arguments - - integer(kind=im), intent(in) :: ncol ! number of layers - integer(kind=im), intent(in) :: nlay ! number of layers - integer(kind=im), intent(in) :: icld ! clear/cloud, cloud overlap flag - integer(kind=im), intent(inout) :: irng ! flag for random number generator - ! 0 = kissvec - ! 1 = Mersenne Twister - integer(kind=im), intent(in) :: nsubcol ! number of sub-columns (g-point intervals) - integer(kind=im), optional, intent(in) :: changeSeed ! allows permuting seed - -! Column state (cloud fraction, cloud water, cloud ice) + variables needed to read physics state - real(kind=rb), intent(in) :: pmid(:,:) ! layer pressure (Pa) - ! Dimensions: (ncol,nlay) -! mji - Add height - real(kind=rb), intent(in) :: hgt(:,:) ! layer height (m) - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: cld(:,:) ! cloud fraction - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: clwp(:,:) ! in-cloud liquid water path (g/m2) - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: ciwp(:,:) ! in-cloud ice water path (g/m2) - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: cswp(:,:) ! in-cloud snow water path (g/m2) - ! Dimensions: (ncol,nlay) - real(kind=rb), intent(in) :: tauc(:,:,:) ! in-cloud optical depth (non-delta scaled) - ! Dimensions: (nbndsw,ncol,nlay) - real(kind=rb), intent(in) :: ssac(:,:,:) ! in-cloud single scattering albedo (non-delta scaled) - ! Dimensions: (nbndsw,ncol,nlay) - real(kind=rb), intent(in) :: asmc(:,:,:) ! in-cloud asymmetry parameter (non-delta scaled) - ! Dimensions: (nbndsw,ncol,nlay) - real(kind=rb), intent(in) :: fsfc(:,:,:) ! in-cloud forward scattering fraction (non-delta scaled) - ! Dimensions: (nbndsw,ncol,nlay) - real(kind=rb), intent(out) :: cld_stoch(:,:,:) ! subcolumn cloud fraction - ! Dimensions: (ngptsw,ncol,nlay) - real(kind=rb), intent(out) :: clwp_stoch(:,:,:) ! subcolumn in-cloud liquid water path - ! Dimensions: (ngptsw,ncol,nlay) - real(kind=rb), intent(out) :: ciwp_stoch(:,:,:) ! subcolumn in-cloud ice water path - ! Dimensions: (ngptsw,ncol,nlay) - real(kind=rb), intent(out) :: cswp_stoch(:,:,:) ! subcolumn in-cloud snow water path - ! Dimensions: (ngptsw,ncol,nlay) - real(kind=rb), intent(out) :: tauc_stoch(:,:,:) ! subcolumn in-cloud optical depth - ! Dimensions: (ngptsw,ncol,nlay) - real(kind=rb), intent(out) :: ssac_stoch(:,:,:) ! subcolumn in-cloud single scattering albedo - ! Dimensions: (ngptsw,ncol,nlay) - real(kind=rb), intent(out) :: asmc_stoch(:,:,:) ! subcolumn in-cloud asymmetry parameter - ! Dimensions: (ngptsw,ncol,nlay) - real(kind=rb), intent(out) :: fsfc_stoch(:,:,:) ! subcolumn in-cloud forward scattering fraction - ! Dimensions: (ngptsw,ncol,nlay) - -! -- Local variables - real(kind=rb) :: cldf(ncol,nlay) ! cloud fraction - ! Dimensions: (ncol,nlay) - -! Mean over the subcolumns (cloud fraction, cloud water , cloud ice) - inactive -! real(kind=rb) :: mean_cld_stoch(ncol,nlay) ! cloud fraction -! real(kind=rb) :: mean_clwp_stoch(ncol,nlay) ! cloud water -! real(kind=rb) :: mean_ciwp_stoch(ncol,nlay) ! cloud ice -! real(kind=rb) :: mean_tauc_stoch(ncol,nlay) ! cloud optical depth -! real(kind=rb) :: mean_ssac_stoch(ncol,nlay) ! cloud single scattering albedo -! real(kind=rb) :: mean_asmc_stoch(ncol,nlay) ! cloud asymmetry parameter -! real(kind=rb) :: mean_fsfc_stoch(ncol,nlay) ! cloud forward scattering fraction - -! Set overlap - integer(kind=im) :: overlap ! 1 = random overlap, 2 = maximum-random, - ! 3 = maximum overlap, 4 = exponential, - ! 5 = exponential-random - real(kind=rb), parameter :: Zo = 2500._rb ! length scale (m) - real(kind=rb), dimension(ncol,nlay) :: alpha ! overlap parameter - -! Constants (min value for cloud fraction and cloud water and ice) - real(kind=rb), parameter :: cldmin = 1.0e-20_rb ! min cloud fraction -! real(kind=rb), parameter :: qmin = 1.0e-10_rb ! min cloud water and cloud ice (not used) - -! Variables related to random number and seed - real(kind=rb), dimension(nsubcol, ncol, nlay) :: CDF, CDF2 ! random numbers - integer(kind=im), dimension(ncol) :: seed1, seed2, seed3, seed4 ! seed to create random number - real(kind=rb), dimension(ncol) :: rand_num ! random number (kissvec) - integer(kind=im) :: iseed ! seed to create random number (Mersenne Twister) - real(kind=rb) :: rand_num_mt ! random number (Mersenne Twister) - -! Flag to identify cloud fraction in subcolumns - logical, dimension(nsubcol, ncol, nlay) :: isCloudy ! flag that says whether a gridbox is cloudy - -! Indices - integer(kind=im) :: ilev, isubcol, i, n, ngbm ! indices - -!------------------------------------------------------------------------------------------ - -! Check that irng is in bounds; if not, set to default - if (irng .ne. 0) irng = 1 - -! Pass input cloud overlap setting to local variable - overlap = icld - -! Ensure that cloud fractions are in bounds - do ilev = 1, nlay - do i = 1, ncol - cldf(i,ilev) = cld(i,ilev) - if (cldf(i,ilev) < cldmin) then - cldf(i,ilev) = 0._rb - endif - enddo - enddo - -! ----- Create seed -------- - -! Advance randum number generator by changeseed values - if (irng.eq.0) then -! For kissvec, create a seed that depends on the state of the columns. Maybe not the best way, but it works. - -! Must use pmid from bottom four layers. - do i=1,ncol - if (pmid(i,1).lt.pmid(i,2)) then - stop 'MCICA_SUBCOL: KISSVEC SEED GENERATOR REQUIRES PMID FROM BOTTOM FOUR LAYERS.' - endif - seed1(i) = (pmid(i,1) - int(pmid(i,1))) * 1000000000_im - seed2(i) = (pmid(i,2) - int(pmid(i,2))) * 1000000000_im - seed3(i) = (pmid(i,3) - int(pmid(i,3))) * 1000000000_im - seed4(i) = (pmid(i,4) - int(pmid(i,4))) * 1000000000_im - enddo - do i=1,changeSeed - call kissvec(seed1, seed2, seed3, seed4, rand_num) - enddo - elseif (irng.eq.1) then - randomNumbers = new_RandomNumberSequence(seed = changeSeed) - endif - - -! ------ Apply overlap assumption -------- - -! generate the random numbers - - select case (overlap) - - - case(1) -! Random overlap -! i) pick a random value at every level - - if (irng.eq.0) then - do isubcol = 1,nsubcol - do ilev = 1,nlay - call kissvec(seed1, seed2, seed3, seed4, rand_num) - CDF(isubcol,:,ilev) = rand_num - enddo - enddo - elseif (irng.eq.1) then - do isubcol = 1, nsubcol - do i = 1, ncol - do ilev = 1, nlay - rand_num_mt = getRandomReal(randomNumbers) - CDF(isubcol,i,ilev) = rand_num_mt - enddo - enddo - enddo - endif - - case(2) -! Maximum-Random overlap -! i) pick a random number for top layer. -! ii) walk down the column: -! - if the layer above is cloudy, we use the same random number than in the layer above -! - if the layer above is clear, we use a new random number - - if (irng.eq.0) then - do isubcol = 1,nsubcol - do ilev = 1,nlay - call kissvec(seed1, seed2, seed3, seed4, rand_num) - CDF(isubcol,:,ilev) = rand_num - enddo - enddo - elseif (irng.eq.1) then - do isubcol = 1, nsubcol - do i = 1, ncol - do ilev = 1, nlay - rand_num_mt = getRandomReal(randomNumbers) - CDF(isubcol,i,ilev) = rand_num_mt - enddo - enddo - enddo - endif - - do ilev = 2,nlay - do i = 1, ncol - do isubcol = 1, nsubcol - if (CDF(isubcol, i, ilev-1) > 1._rb - cldf(i,ilev-1) ) then - CDF(isubcol,i,ilev) = CDF(isubcol,i,ilev-1) - else - CDF(isubcol,i,ilev) = CDF(isubcol,i,ilev) * (1._rb - cldf(i,ilev-1)) - endif - enddo - enddo - enddo - - - case(3) -! Maximum overlap -! i) pick same random numebr at every level - - if (irng.eq.0) then - do isubcol = 1,nsubcol - call kissvec(seed1, seed2, seed3, seed4, rand_num) - do ilev = 1,nlay - CDF(isubcol,:,ilev) = rand_num - enddo - enddo - elseif (irng.eq.1) then - do isubcol = 1, nsubcol - do i = 1, ncol - rand_num_mt = getRandomReal(randomNumbers) - do ilev = 1, nlay - CDF(isubcol,i,ilev) = rand_num_mt - enddo - enddo - enddo - endif - - -! mji - Activate exponential cloud overlap option - case(4) - ! Exponential overlap: weighting between maximum and random overlap increases with the distance. - ! The random numbers for exponential overlap verify: - ! j=1 RAN(j)=RND1 - ! j>1 if RND1 < alpha(j,j-1) => RAN(j) = RAN(j-1) - ! RAN(j) = RND2 - ! alpha is obtained from the equation - ! alpha = exp(-(Z(j)-Z(j-1))/Zo) where Zo is a characteristic length scale - - ! compute alpha - do i = 1, ncol - alpha(i, 1) = 0._rb - do ilev = 2,nlay - alpha(i, ilev) = exp( -( hgt (i, ilev) - hgt (i, ilev-1)) / Zo) - enddo - enddo - - ! generate 2 streams of random numbers - if (irng.eq.0) then - do isubcol = 1,nsubcol - do ilev = 1,nlay - call kissvec(seed1, seed2, seed3, seed4, rand_num) - CDF(isubcol, :, ilev) = rand_num - call kissvec(seed1, seed2, seed3, seed4, rand_num) - CDF2(isubcol, :, ilev) = rand_num - enddo - enddo - elseif (irng.eq.1) then - do isubcol = 1, nsubcol - do i = 1, ncol - do ilev = 1, nlay - rand_num_mt = getRandomReal(randomNumbers) - CDF(isubcol,i,ilev) = rand_num_mt - rand_num_mt = getRandomReal(randomNumbers) - CDF2(isubcol,i,ilev) = rand_num_mt - enddo - enddo - enddo - endif - - ! generate random numbers - do ilev = 2,nlay - where (CDF2(:, :, ilev) < spread(alpha (:,ilev), dim=1, nCopies=nsubcol) ) - CDF(:,:,ilev) = CDF(:,:,ilev-1) - end where - end do - -! mji - Activate exponential-random cloud overlap option - case(5) - ! Exponential-random overlap: -! call wrf_error_fatal("Cloud Overlap case 5: ER has not yet been implemented. Stopping...") - - end select - - -! -- generate subcolumns for homogeneous clouds ----- - do ilev = 1, nlay - isCloudy(:,:,ilev) = (CDF(:,:,ilev) >= 1._rb - spread(cldf(:,ilev), dim=1, nCopies=nsubcol) ) - enddo - -! where the subcolumn is cloudy, the subcolumn cloud fraction is 1; -! where the subcolumn is not cloudy, the subcolumn cloud fraction is 0; -! where there is a cloud, define the subcolumn cloud properties, -! otherwise set these to zero - - ngbm = ngb(1) - 1 - do ilev = 1,nlay - do i = 1, ncol - do isubcol = 1, nsubcol - if ( iscloudy(isubcol,i,ilev) ) then - cld_stoch(isubcol,i,ilev) = 1._rb - clwp_stoch(isubcol,i,ilev) = clwp(i,ilev) - ciwp_stoch(isubcol,i,ilev) = ciwp(i,ilev) - cswp_stoch(isubcol,i,ilev) = cswp(i,ilev) - n = ngb(isubcol) - ngbm - tauc_stoch(isubcol,i,ilev) = tauc(n,i,ilev) - ssac_stoch(isubcol,i,ilev) = ssac(n,i,ilev) - asmc_stoch(isubcol,i,ilev) = asmc(n,i,ilev) - fsfc_stoch(isubcol,i,ilev) = fsfc(n,i,ilev) - else - cld_stoch(isubcol,i,ilev) = 0._rb - clwp_stoch(isubcol,i,ilev) = 0._rb - ciwp_stoch(isubcol,i,ilev) = 0._rb - cswp_stoch(isubcol,i,ilev) = 0._rb - tauc_stoch(isubcol,i,ilev) = 0._rb - ssac_stoch(isubcol,i,ilev) = 1._rb - asmc_stoch(isubcol,i,ilev) = 0._rb - fsfc_stoch(isubcol,i,ilev) = 0._rb - endif - enddo - enddo - enddo - - -! -- compute the means of the subcolumns --- -! mean_cld_stoch(:,:) = 0._rb -! mean_clwp_stoch(:,:) = 0._rb -! mean_ciwp_stoch(:,:) = 0._rb -! mean_tauc_stoch(:,:) = 0._rb -! mean_ssac_stoch(:,:) = 0._rb -! mean_asmc_stoch(:,:) = 0._rb -! mean_fsfc_stoch(:,:) = 0._rb -! do i = 1, nsubcol -! mean_cld_stoch(:,:) = cld_stoch(i,:,:) + mean_cld_stoch(:,:) -! mean_clwp_stoch(:,:) = clwp_stoch( i,:,:) + mean_clwp_stoch(:,:) -! mean_ciwp_stoch(:,:) = ciwp_stoch( i,:,:) + mean_ciwp_stoch(:,:) -! mean_tauc_stoch(:,:) = tauc_stoch( i,:,:) + mean_tauc_stoch(:,:) -! mean_ssac_stoch(:,:) = ssac_stoch( i,:,:) + mean_ssac_stoch(:,:) -! mean_asmc_stoch(:,:) = asmc_stoch( i,:,:) + mean_asmc_stoch(:,:) -! mean_fsfc_stoch(:,:) = fsfc_stoch( i,:,:) + mean_fsfc_stoch(:,:) -! end do -! mean_cld_stoch(:,:) = mean_cld_stoch(:,:) / nsubcol -! mean_clwp_stoch(:,:) = mean_clwp_stoch(:,:) / nsubcol -! mean_ciwp_stoch(:,:) = mean_ciwp_stoch(:,:) / nsubcol -! mean_tauc_stoch(:,:) = mean_tauc_stoch(:,:) / nsubcol -! mean_ssac_stoch(:,:) = mean_ssac_stoch(:,:) / nsubcol -! mean_asmc_stoch(:,:) = mean_asmc_stoch(:,:) / nsubcol -! mean_fsfc_stoch(:,:) = mean_fsfc_stoch(:,:) / nsubcol - - end subroutine generate_stochastic_clouds_sw - - -!-------------------------------------------------------------------------------------------------- - subroutine kissvec(seed1,seed2,seed3,seed4,ran_arr) -!-------------------------------------------------------------------------------------------------- - -! public domain code made available from http://www.fortran.com/ -! downloaded by pjr on 03/16/04 for NCAR CAM -! converted to vector form, functions inlined by pjr,mvr on 05/10/2004 - -! The KISS (Keep It Simple Stupid) random number generator. Combines: -! (1) The congruential generator x(n)=69069*x(n-1)+1327217885, period 2^32. -! (2) A 3-shift shift-register generator, period 2^32-1, -! (3) Two 16-bit multiply-with-carry generators, period 597273182964842497>2^59 -! Overall period>2^123; - -! - real(kind=rb), dimension(:), intent(inout) :: ran_arr - integer(kind=im), dimension(:), intent(inout) :: seed1,seed2,seed3,seed4 - integer(kind=im) :: i,sz,kiss - integer(kind=im) :: m, k, n - -! inline function - m(k, n) = ieor (k, ishft (k, n) ) - - sz = size(ran_arr) - do i = 1, sz - seed1(i) = 69069_im * seed1(i) + 1327217885_im - seed2(i) = m (m (m (seed2(i), 13_im), - 17_im), 5_im) - seed3(i) = 18000_im * iand (seed3(i), 65535_im) + ishft (seed3(i), - 16_im) - seed4(i) = 30903_im * iand (seed4(i), 65535_im) + ishft (seed4(i), - 16_im) - kiss = seed1(i) + seed2(i) + ishft (seed3(i), 16_im) + seed4(i) - ran_arr(i) = kiss*2.328306e-10_rb + 0.5_rb - end do - - end subroutine kissvec - -!! @} - -! !........................................! end module rrtmg_sw ! !========================================! diff --git a/physics/radsw_main.meta b/physics/radsw_main.meta index 187d26f21..632dc3912 100644 --- a/physics/radsw_main.meta +++ b/physics/radsw_main.meta @@ -1,7 +1,7 @@ [ccpp-table-properties] name = rrtmg_sw type = scheme - dependencies = mersenne_twister.f,physcons.F90,physparam.f,radsw_datatb.f,radsw_param.f,HWRF_mcica_random_numbers.F90,HWRF_mersenne_twister.F90 + dependencies = mersenne_twister.f,physcons.F90,physparam.f,radsw_datatb.f,radsw_param.f ######################################################################## [ccpp-arg-table]