!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! GNU General Public License !! !! !! !! This file is part of the Flexible Modeling System (FMS). !! !! !! !! FMS is free software; you can redistribute it and/or modify !! !! it and are expected to follow the terms of the GNU General Public !! !! License as published by the Free Software Foundation. !! !! !! !! FMS is distributed in the hope that it will be useful, !! !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! !! GNU General Public License for more details. !! !! !! !! You should have received a copy of the GNU General Public License !! !! along with FMS; if not, write to: !! !! Free Software Foundation, Inc. !! !! 59 Temple Place, Suite 330 !! !! Boston, MA 02111-1307 USA !! !! or see: !! !! http://www.gnu.org/licenses/gpl.txt !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! MODULE DIAG_CLOUD_MOD !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! DIAGNOSTIC CLOUD PREDICTION - Gordon (1992) ! ! 1999 Feb -> 2000 July ! Contact persons: Bill Stern (for code structure information) ! Tony Gordon (for cloud scheme information) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !------------------------------------------------------------------- ! Calculates cloud fractions diagnostically using relative humidity, ! omega and stability !------------------------------------------------------------------- use fms_mod, only: error_mesg, FATAL, NOTE, file_exist, & check_nml_error, open_namelist_file, & mpp_pe, mpp_root_pe, close_file, & read_data, write_data, & write_version_number, stdlog, open_restart_file use fms_io_mod, only: register_restart_field, restart_file_type, & save_restart, restore_state, get_restart_io_mode use Constants_Mod, only: Cp_Air, rdgas, rvgas, Kappa, HLv use time_manager_mod, only: TIME_TYPE use cloud_zonal_mod, only: CLOUD_ZONAL_INIT, GETCLD use diag_cloud_rad_mod, only: CLOUD_TAU_DRIVER, diag_cloud_rad_INIT,& cloud_pres_thick_for_tau, & cloud_opt_prop_tg_lw, & cloud_opt_prop_tg_sw, & cloud_optical_depths, & cloud_optical_depths2 use sat_vapor_pres_mod, ONLY: compute_qs use shallow_conv_mod, ONLY: SHALLOW_CONV_INIT,MYLCL !----------------------------------------------------------------------- implicit none !----------------------------------------------------------------------- private !--------------------- version number ---------------------------------- character(len=128) :: version = '$Id: diag_cloud.F90,v 17.0 2009/07/21 02:54:07 fms Exp $' character(len=128) :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $' logical :: module_is_initialized = .false. !----------------------------------------------------------------------- ! parmameter mxband = max number of radiative bands to be considered for some ! cloud properties integer, parameter :: mxband = 4 ! parmameter nrhc = number of critical relative humidity threshold values ! as a function of vertical layers integer, parameter :: nrhc = 3 ! **************************************************************************** ! WARNING: The parameter mxband and the namelist variable nband should have ! the same value. Is there a way to guarantee this or can an error check ! be added to the code to make sure that this condition is NOT violated? ! **************************************************************************** ! The following constants are derive from other constants in constants_mod real, parameter :: d622 = rdgas/rvgas, d378 = 1.0-d622 real, parameter :: p00 = 1000.0E2 !----------------------------------------------------------------------- ! **** parameter used and defined in def_hi_mid_low **** ! real, parameter :: trshld_camt = 0.25 ! TRSHLD_CAMT - This is a cloud amounht threshold value, used in ! conjunction with the thick cloud namelist options. If a thick ! cloud option is on then the level above the cloud amount max ! for a particular cloud type is considered to be part of an ! extended cloud if it is within this fraction of the max cloud amount. ! !----------------- arrays for cloud predictor averaging code -------------------------- real, allocatable, dimension (:,:,:) :: temp_sum,qmix_sum,rhum_sum real, allocatable, dimension (:,:) :: qmix_sum2 real, allocatable, dimension (:,:,:) :: omega_sum real, allocatable, dimension (:,:,:) :: lgscldelq_sum,cnvcntq_sum real, allocatable, dimension (:,:) :: convprc_sum integer, allocatable, dimension (:,:) :: nsum, nsum2 !----------------------------------------------------------------------- ! for netcdf restart type(restart_file_type), save :: Dia_restart !--------------------------------------------------------------------- ! logical switches for initialization calls: ! if = .false. -> has not been called logical :: do_cpred_init = .false. logical :: do_crad_init = .false. !--------------------------------------------------------------------- ! --- NAMELIST (diag_cloud_nml) !--------------------------------------------------------------------- ! RHC - critical humidity value (ras = 0.8 - 0.84, mca =0.7) ! (note: in vers >= 0.9.3 a function of 3 levels ! i.e.,"high" , "mid", "low" - but here is more general) ! PBOUNDS - sets pressure bounds for RHC (dimension = size(rhc) - 1 ! DO_AVERAGE - logical flag for time averaging cloud predictor variables ! LQUADRA - logical switch for turning on quadratic relation ! for calculating rhum clouds from rhum, ! i.e., true for quadratice scheme, false for linear scheme ! LRHCNV - logical switch for using rhum fields as follows: ! if true - use rel humidities modified for presence of ! convective clouds (rhumcnv), otherwise use original ! rel humidities (rhum) ! LOMEGA - logical switch for turning on omega correction to rhum ! clouds - true for omega correction, otherwise false ! LCNVCLD - logical switch for turning on calculation of deep convective ! clouds - true for deep convective clouds, otherwise false ! L_THEQV - logical switch for turning on calculation of shallow convective ! clouds - true for shallow convective clouds, otherwise false ! LINVERS - logical switch for turning on calculation of marine stratus ! clouds - true for marine stratus, otherwise false ! LSLINGO - logical variable = true apply Slingo marine stratus ! scheme, otherwise = false. ! LREGRSC - logical variable = true apply Tim Li marine stratus ! scheme, otherwise = false. Slingo & Li schemes may be ! used in combination, but atleast one scheme must be used. ! LTHICK_HIGH - logical variable = true -> allow possibility of raising ! high cloud tops one sigma level to increase their thickness ! from 1 to nmax levels; otherwise they remain thin ! (1 level) ! LTHICK_MID - logical variable = true -> allow possibility of raising ! mid cloud tops one sigma level to increase their thickness ! from 1 to nmax levels; otherwise they remain thin ! (1 level) ! LTHICK_LOW - logical variable = true -> allow possibility of raising ! low cloud tops one sigma level to increase their thickness ! from 1 to nmax levels; otherwise they remain thin ! (1 level) ! HIGH_LEV_CLOUD_INDEX - level above which no clouds are allowed to form ! (model level index) ! NOFOG - logical switch for not allowing rhum clouds (or fog) ! to occur beneath a certain level (low_lev_cloud_index) -> ! nofog = true ! to allow clouds at the lowest model level -> nofog = false ! LOW_LEV_CLOUD_INDEX - level below which no clouds are allowed to occur ! when nofog = true (model level index) ! NBAND - max number of radiative bands to be considered for some ! cloud properties ! PSHALLOW - top of shallow convective layer (pressure level - n/m**2 ) ! WCUT0 - omega cutoff value for omega cloud depletion factor = 0 ! WCUT1 - omega cutoff value for omega cloud depletion factor = 1 ! t_cold temperature defining ice-liquid cloud transition !----------------------------------------------------------------------- real :: t_cold= 263.16 real :: & pshallow=.750E+05, wcut0 = .10, wcut1 = 0.0 real, dimension(nrhc) :: rhc = (/ 0.8,0.8,0.84 /) real, dimension(nrhc-1) :: pbounds = (/ .400E5, .750e5 /) integer :: & high_lev_cloud_index=3, low_lev_cloud_index=16, nband=4 logical :: & do_average = .true., & lquadra = .true., nofog = .false.,lrhcnv = .false., & lomega = .true.,lcnvcld = .true.,l_theqv = .true., & linvers = .false.,lslingo = .true., lregrsc = .true., & lthick_high = .true.,lthick_mid = .true.,lthick_low = .true. NAMELIST / diag_cloud_nml / & rhc,pbounds,do_average,lquadra,lrhcnv,lomega,lcnvcld,l_theqv, & linvers,lslingo,lregrsc,lthick_high,lthick_mid,lthick_low, & high_lev_cloud_index, nofog, low_lev_cloud_index, nband, & pshallow, wcut0, wcut1, t_cold integer :: num_pts, tot_pts logical :: do_netcdf_restart public diag_cloud_driver, diag_cloud_init, diag_cloud_end public diag_cloud_driver2 public diag_cloud_sum, diag_cloud_avg, diag_cloud_avg2, do_diag_cloud public diag_cloud_restart contains !############################################################################# SUBROUTINE DIAG_CLOUD_DRIVER (is,js, & temp,qmix,rhum,omega,lgscldelq,cnvcntq,convprc, & pfull,phalf,psfc,coszen,lat,time, & nclds,cldtop,cldbas,cldamt,r_uv,r_nir,ab_uv,ab_nir, & em_lw, conc_drop, conc_ice, size_drop, size_ice, & kbot) ! Arguments (intent in) integer, intent(in) :: is,js type(time_type), intent(in) :: time real, intent(in) :: lat(:,:) real, intent(in), dimension (:,:,:) :: temp,qmix,rhum,omega real, intent(in), dimension (:,:,:) :: lgscldelq,cnvcntq,pfull, phalf real, intent(in), dimension (:,:) :: convprc,psfc, coszen integer, intent(in), OPTIONAL, dimension(:,:) :: kbot ! INPUT ! ------ ! IS,JS starting i,j indices from the full horizontal grid ! IX, IY Horizontal dimensions for global storage arrays ! TEMP Temperature (Deg K) at full model levels ! (dimensioned IDIM x JDIM x kx) ! QMIX Mixing Ratio at full model levels ! (dimensioned IDIM x JDIM x kx) ! RHUM Relative humidity fraction at full model levels ! (dimensioned IDIM x JDIM x kx) ! OMEGA Pressure vertical velocity at full model levels ! (dimensioned IDIM x JDIM x kx) ! LGSCLDELQ Averaged rate of change in mix ratio due to lg scale precip ! at full model levels ! (dimensioned IDIM x JDIM x kx) ! CNVCNTQ Accumulated count of change in mix ratio due to conv precip ! at full model levels ! (dimensioned IDIM x JDIM x kx) ! convprc Accumulated conv precip rate summed over all ! full model levels (mm/day ) ! (dimensioned IDIM x JDIM) ! PFULL Pressure at full model levels ! (dimensioned IDIM x JDIM x kx) ! PHALF Pressure at half model levels ! (dimensioned IDIM x JDIM x kx+1) ! PSFC Surface pressure field ! (dimensioned IDIM x JDIM) ! COSZEN cosine of the zenith angle ! (dimensioned IDIM x JDIM) ! TIME time of year (time_type) ! LAT latitudes in radians, dimensioned by (1xJDIM) ! KBOT OPTIONAL; lowest model level index array ! (dimensioned IDIM x JDIM) !=================================================================== ! Arguments (intent out) integer, intent(out), dimension(:,:,:) :: cldtop,cldbas integer, intent(out), dimension(:,:) :: nclds real, intent(out), dimension(:,:,:), optional :: r_uv,r_nir,ab_uv, & ab_nir,em_lw, & conc_drop, conc_ice, & size_drop, size_ice real, intent(out), dimension(:,:,:) :: cldamt ! OUTPUT ! ------ ! NCLDS number of (random overlapping) clouds in column and also ! the current # for clouds to be operating on ! (dimensioned IDIM x JDIM ) ! CLDTOP index of cloud tops (at cloud levels) ! (dimensioned IDIM x JDIM x kx) ! CLDBAS index of cloud bottoms (at cloud levels) ! (dimensioned IDIM x JDIM x kx) ! CLDAMT cloud amount (fraction) (at cloud levels) ! (dimensioned IDIM x JDIM x kx) ! R_UV fractional amount of ultraviolet radiation ! reflected by the clouds (at cloud levels) ! R_NIRfractional amount of near inrared radiation ! reflected by the clouds (at cloud levels) ! AB_UVfractional amount of ultraviolet radiation ! absorbed by the clouds (at cloud levels) ! AB_NIRfractional amount of near inrared radiation ! absorbed by the clouds (at cloud levels) ! EM_LWemissivity for the clouds (at cloud levels) !======================================================================= ! (Intent local) integer, dimension(size(rhum,1),size(rhum,2)) :: & lhight, lhighb, lmidt, lmidb, llowt integer, dimension(size(rhum,1),size(rhum,2),size(rhum,3)) :: icld real, dimension(size(rhum,1),size(rhum,2)) :: qmix_kx real, dimension(size(rhum,1),size(rhum,2),size(rhum,3),mxband) :: tau real, dimension(size(rhum,1),size(rhum,2),size(rhum,3)) :: & tempcld,delp_true, delp integer idim, jdim, kx, lk logical :: rad_prop, wat_prop integer :: max_cld integer :: i,j,k,n ! TAU cloud optical depth (at cloud levels) ! (dimensioned IDIM x JDIM x kx x MXBAND) ! TEMPCLD cloud layer mean temperature (degrees Kelvin) ! (at cloud levels) ! DELP_TRUE true cloud pressure thickness of distinct cloud layers ! (at cloud levels) ! QMIX_KX Lowest level mixing ratio ! (dimensioned IDIM x JDIM) ! ICLD marker array of cloud types/heights (at cloud levels) ! (dimensioned IDIM x JDIM x kx) ! LHIGHT vertical level index upper limit for high cloud tops ! a function of lon and latitude (dimensioned IDIMxJDIM) ! LHIGHB vertical level index lower limit for high cloud bases ! a function of lon and latitude (dimensioned IDIMxJDIM) ! LMIDT vertical level index upper limit for mid cloud tops ! a function of lon and latitude (dimensioned IDIMxJDIM) ! LMIDB vertical level index lower limit for mid cloud bases ! a function of lon and latitude (dimensioned IDIMxJDIM) ! LLOWT vertical level index upper limit for low cloud tops ! a function of lon and latitude (dimensioned IDIMxJDIM) ! IERR Error flag ! IPRNT Longitude index for sample printing ! ! note: LK - vertical level below which no low cloud bases can exist, is ! calculated in routines where needed using namelist inputs ! low_lev_cloud_index and nofog for the cloud amount prediction code ! but is passed to the cloud raidiative properties code as an argument. !=================================================================== !----------------------------------------------------------------------- idim = size(rhum,1) jdim = size(rhum,2) kx = size(rhum,3) tau = 0.0 rad_prop = .false. wat_prop = .false. if (present (r_uv) .or. present(r_nir) .or. present (ab_uv) & .or. present(ab_nir) .or. present (em_lw) ) then rad_prop = .true. r_uv = 0. r_nir = 0. ab_uv = 0. ab_nir = 0. em_lw = 0. endif if (present(conc_drop) .or. present(conc_ice) .or. & present(size_drop) .or. present(size_ice) ) then wat_prop = .true. endif if ( (.not. rad_prop) .and. (.not. wat_prop) ) then rad_prop = .true. endif ! define lower limit for low cloud bases if (nofog) then lk = low_lev_cloud_index else lk = kx endif ! cldtim drives cloud prediction scheme call cldtim ( temp,qmix,rhum,omega,lgscldelq,cnvcntq,convprc, & pfull, phalf,psfc, lat, time, tempcld,delp_true, & cldtop,cldbas,cldamt, & lhight,lhighb, lmidt, lmidb, llowt,icld,nclds,kbot) ! lowest level mixing ratio for anomalous absorption in cloud-radiation qmix_kx(:,:) = qmix(:,:,kx) max_cld = MAXVAL(nclds(:,:)) IF (max_cld .gt. 0) then call cloud_pres_thick_for_tau (nclds,icld,cldtop,cldbas, & & delp_true,lhight,lhighb, lmidt, lmidb, llowt,lk, delp, & & phalf, psfc ) call cloud_optical_depths(nclds,icld,cldtop,cldbas,tempcld,delp, & tau,phalf ) ! cloud_tau_driver drives cloud radiative properties scheme if (rad_prop) then call cloud_tau_driver ( qmix_kx, & tempcld, & tau, coszen, & r_uv=r_uv, r_nir=r_nir, ab_nir=ab_nir, ab_uv=ab_uv, & em_lw=em_lw) endif endif ! (max_cld > 0) !----------------------------------------------------------------------- ! print output for 1-d testing ! iprnt = 1 ! print *, ' diag_cloud sample output for iprnt,js= ', iprnt,js ! print *,'cloud layers = ', nclds(iprnt,1) ! print *,'cloud amounts = ', cldamt(iprnt,1,:) ! print *,'cloud tops = ', cldtop(iprnt,1,:) ! print *,'cloud bases = ', cldbas(iprnt,1,:) ! print *,'cloud types = ', icld(iprnt,1,:) ! print *, ' clouds_rad_tg sample output ' ! print *,'reflectivity uv = ', r_uv(iprnt,1,:) ! print *,'reflectivity nir = ', r_nir(iprnt,1,:) ! print *,'absorptivity uv = ', ab_uv(iprnt,1,:) ! print *,'absorptivity nir = ', ab_nir(iprnt,1,:) ! print *,'emissivity = ', em_lw(iprnt,1,:) ! print *,'optical depth = ', tau(iprnt,1,:,:) !----------------------------------------------------------------------- if (present (conc_drop)) then !----------------------------------------------------------------------- !!!RSH !! NOTE: ! THE FOLLOWING is here as an INITIAL IMPLEMENTATION to allow compil- ! ation and model execution, and provide "reasonable ?? " values. ! Code developed but NOT YET ADDED HERE reflects the current approach. ! That code is available under the fez release, and will be added to ! the repository when upgrades to the cloud-radiation modules are com- ! pleted. !!!RSH ! obtain drop and ice size and concentration here, consistent with the ! diag_cloud scheme. ! As a test case, ! the following is a simple specification of constant concentration and ! size in all boxes defined as cloudy, attempting to come close to ! the prescribd values in microphys_rad. ! assume ice cld thickness = 2.0 km; then conc_ice=10.0E-03 => ! iwp = 20 g/m^2, similar to that prescribed in microphys_rad. ! assume water cld thickness = 3.5 km; then conc_drop = 20E-03 => ! lwp = 70 g / m^2, similar to that prescribed in microphys_rad. ! use sizes as used in microphys_rad (50 and 20 microns). when done, ! radiative boundary fluxes are "similar" to non-microphysical results ! for test case done here, and shows reasonable sensitivity to ! variations in concentrations. conc_ice = 0. conc_drop = 0. size_ice = 50. size_drop = 20. IF (max_cld .gt. 0) then idim = size(tau,1) jdim = size(tau,2) do j= 1,jdim do i=1,idim do n=1,nclds(i,j) do k=cldtop(i,j,n), cldbas(i,j,n) if (tempcld(i,j,n) < t_cold) then conc_ice(i,j,k) = 10.0E-03 ! units : g/m^3 size_ice(i,j,k) = 50. ! units : diameter in microns else conc_drop(i,j,k) = 20.0E-03 ! units : g/m^3 size_drop(i,j,k) = 20. ! units : diameter in microns endif end do end do end do end do endif endif end SUBROUTINE DIAG_CLOUD_DRIVER !--------------------------------------------------------------------- subroutine diag_cloud_driver2 (is, js, press, pflux, lat, time, nclds, & cldtop, cldbas, cldamt, liq_frac, tau, & ice_cloud, kbot) !--------------------------------------------------------------------- ! diag_cloud_driver2 returns the cloud specification arrays for the ! gordon diag cloud scheme. returned are the number of clouds per ! column, the cloud top, cloud base and fractional coverage of each ! cloud, the amount of the cloud which is liquid, its optical depth ! and an indicator as to whether it is ice or liquid (a different ! criterion than is used for the liquid fraction determination). !---------------------------------------------------------------------- integer, intent(in) :: is,js real, dimension (:,:,:), intent(in) :: press, pflux real, dimension(:,:), intent(in) :: lat type(time_type), intent(in) :: time integer, dimension(:,:), intent(inout) :: nclds integer, dimension(:,:,:), intent(out) :: cldtop,cldbas real, dimension(:,:,:), intent(out) :: cldamt, liq_frac real, dimension(:,:,:,:), intent(out) :: tau logical, dimension(:,:,:), intent(out) :: ice_cloud integer, dimension(:,:), intent(in), optional :: kbot !--------------------------------------------------------------------- ! intent(in) variables: ! ! is,js starting subdomain i,j indices of data ! in the physics_window being integrated ! press pressure at model levels (1:nlev), surface ! pressure is stored at index value nlev+1 ! [ (kg /( m s^2) ] ! pflux average of pressure at adjacent model levels ! [ (kg /( m s^2) ] ! lat latitude of model points [ radians ] ! time time at which radiation calculation is to apply ! [ time_type (days, seconds) ] ! ! intent(inout) variables: ! ! nclds total number of clouds in each grid column ! ! intent(out) variables: ! ! cldtop k index of cloud top for each cloud ! cldbas k index of cloud base for each cloud ! cldamt fractional cloudiness for each cloud ! [ dimensionless ] ! liq_frac fraction of cloud which is liquid ! [ dimensionless ] ! tau cloud optical depth [ dimensionless ] ! ice_cloud logical flag indicating whether cloud is liquid ! or ice ! ! intent(in), optional variables: ! ! kbot present when running eta vertical coordinate, ! index of lowest model level above ground ! !--------------------------------------------------------------------- !-------------------------------------------------------------------- ! local variables: integer, dimension (size(press,1), size(press,2)) :: & lhight, lhighb, lmidt, lmidb, llowt integer, dimension (size(press,1), size(press,2), & size(press,3)-1) :: icld real, dimension (size(press,1), size(press,2), & size(press,3)-1) :: & temp, qmix, rhum, omega, lgscldelq, cnvcntq, & delp, tempcld, delp_true, pfull real, dimension (size(press,1), size(press,2)) :: convprc, psfc real, dimension(size(press,1), size(press,2), & size(press,3)) :: phalf integer :: kx, lk, ierr, max_cld integer :: i, j, n !--------------------------------------------------------------------- ! local variables: ! ! lhight vertical level index upper limit for high cloud tops ! a function of lat and lon ! lhighb vertical level index lower limit for high cloud bases ! a function of lat and lon ! lmidt vertical level index upper limit for mid cloud tops ! a function of lat and lon ! lmidb vertical level index lower limit for mid cloud bases ! a function of lat and lon ! llowt vertical level index upper limit for low cloud tops ! a function of ltt and lon ! icld marker array of cloud types/heights (at cloud levels) ! temp temperature at full model levels [ deg K ] ! qmix mixing ratio at full model levels [ kg H2O / kg air ] ! rhum relative humidity fraction at full model levels ! [ dimensionless ] ! omega pressure vertical velocity at full model levels ! [ mb / sec ??????? ] ! lgscldelq averaged rate of change in mixing ratio due to large ! scale precip at full model levels ! cnvcntq accumulated count of change in mixing ratio due to ! convective precip at full model levels ! delp pressure thickness of model layers [ kg / (m s^2) ] ! tempcld cloud layer mean temperature, at cloud levels ! [ deg K ] ! delp_true true cloud pressure thickness of distinct cloud ! layers (at cloud levels) [ kg / (m s^2) ] ! pfull pressure at full levels [ kg / (m s^2) ] ! convprc accumulated conv precip rate summed over all ! full model levels [ mm/day ] ! psfc surface pressure field [ kg / (m s^2) ] ! phalf pressure at model half levels [ kg / (m s^2) ] ! kx number of model layers ! lk vertical level below which no low cloud bases can ! exist ! ierr error flag ! max_cld max number of clouds in any column in the current ! physics window ! i,j,n do loop indices ! !--------------------------------------------------------------------- !---------------------------------------------------------------------- ! define the number of model layers. !---------------------------------------------------------------------- kx = size(press,3) - 1 !--------------------------------------------------------------------- ! define the needed pressure arrays. !--------------------------------------------------------------------- pfull(:,:,:) = press(:,:,1:kx) phalf(:,:,:) = pflux(:,:,:) psfc(:,:) = press(:,:,kx+1) !-------------------------------------------------------------------- ! call diag_cloud_avg to obtain the appropriate values for the input ! arrays needed to define the cloud locations and amounts. these may ! or may not be time-averaged values. !--------------------------------------------------------------------- call diag_cloud_avg (is, js, temp, qmix, rhum, omega, lgscldelq, & cnvcntq, convprc, ierr) !---------------------------------------------------------------------- ! initialize the output fields produced by this module. !--------------------------------------------------------------------- tau = 0. liq_frac = 0. cldamt = 0. cldtop = 0 cldbas = 0 ice_cloud = .false. !---------------------------------------------------------------------- ! if input data was appropriately returned from diag_cloud_avg, ! proceed with the determination of the cloud field. !--------------------------------------------------------------------- if (ierr == 0) then !--------------------------------------------------------------------- ! define the lowest model level which can be a cloud base. it is ! either the lowest model level, or a level determined from namelist ! input. !--------------------------------------------------------------------- if (nofog) then lk = low_lev_cloud_index else lk = kx endif !-------------------------------------------------------------------- ! call cldtim to drive the cloud prediction scheme. !-------------------------------------------------------------------- call cldtim (temp, qmix, rhum, omega, lgscldelq, cnvcntq, & convprc, pfull, phalf, psfc, lat, time, tempcld, & delp_true, cldtop, cldbas, cldamt, lhight, lhighb,& lmidt, lmidb, llowt, icld, nclds, kbot) !--------------------------------------------------------------------- ! determine the maximum number of clouds in any of the columns in the ! physics window. !--------------------------------------------------------------------- max_cld = MAXVAL(nclds(:,:)) !-------------------------------------------------------------------- ! if cloud is present anywhere in the window, call ! cloud_pres_thick_for_tau to determine the cloud thicknesses and ! the call cloud_optical_depths2 to determine the optical depths and ! liquid fraction of each cloud. !--------------------------------------------------------------------- if (max_cld > 0) then call cloud_pres_thick_for_tau (nclds, icld, cldtop, cldbas, & delp_true, lhight, lhighb, & lmidt, lmidb, llowt, lk, delp,& phalf, psfc) call cloud_optical_depths2 (nclds, icld, cldtop, cldbas, & tempcld, delp, tau, phalf, & liq_frac) endif !--------------------------------------------------------------------- ! determine whether the cloud temperature will support a liquid or ! an ice cloud. the parameter t_cold is the cutoff temperature value. !--------------------------------------------------------------------- do j= 1,size(press,2) do i=1,size(press,1) do n=1,nclds(i,j) if (tempcld(i,j,n) < t_cold) then ice_cloud(i,j,n) = .true. endif end do end do end do !---------------------------------------------------------------------- ! if input data was not acceptably returned from diag_cloud_avg, ! determine if this represents an error, or is just a manifestation ! of coldstart behavior. if this is coldstart step, set clouds to ! zero and continue. !--------------------------------------------------------------------- else if (num_pts >= tot_pts) then call error_mesg ('diag_cloud_mod', & ' no diag cloud data available; ierr /= 0', FATAL) else num_pts = num_pts + size(press,1)*size(press,2) nclds = 0 endif endif ! (ierr=0) !-------------------------------------------------------------------- end subroutine diag_cloud_driver2 !--------------------------------------------------------------------- !################################################################## SUBROUTINE CLDTIM (temp,qmix,rhum, omega,lgscldelq,cnvcntq,convprc, & pfull, phalf,psfc, lat, time, tempcld,delp_true, & cldtop,cldbas,cldamt, & lhight,lhighb, lmidt, lmidb, llowt,icld,nclds,kbot) !=================================================================== ! Arguments (intent in) ! Namelist variables used in this routine (defined at top of module) ! HIGH_LEV_CLOUD_INDEX - level above which no clouds are allowed to form ! (model level index) ! NOFOG - logical switch for not allowing rhum clouds (or fog) ! to occur beneath a certain level (low_lev_cloud_index) -> ! nofog = true ! to allow clouds at the lowest model level -> nofog = false ! LOW_LEV_CLOUD_INDEX - level below which no clouds are allowed to occur ! when nofog = true (model level index) type(time_type), intent(in) :: time real, intent(in) :: lat(:,:) real, intent(in), dimension (:,:,:) :: temp,qmix,rhum,omega real, intent(in), dimension (:,:,:) :: lgscldelq,cnvcntq,pfull, phalf real, intent(in), dimension (:,:) :: convprc,psfc ! ! INPUT ! ----- ! ! TEMP Temperature (Deg K) at full model levels ! (dimensioned IDIM x JDIM x kx) ! QMIX Mixing Ratio at full model levels ! (dimensioned IDIM x JDIM x kx) ! RHUM Relative humidity fraction at full model levels ! (dimensioned IDIM x JDIM x kx) ! OMEGA Pressure vertical velocity at full model levels ! (dimensioned IDIM x JDIM x kx) ! LGSCLDELQ Averaged rate of change in mix ratio due to lg scale precip ! at full model levels ! (dimensioned IDIM x JDIM x kx) ! CNVCNTQ Accumulated count of change in mix ratio due to conv precip ! at full model levels ! (dimensioned IDIM x JDIM x kx) ! CONVPRC Accumulated conv precip rate summed over all ! full model levels (mm/day ) ! (dimensioned IDIM x JDIM) ! PFULL Pressure at full model levels ! (dimensioned IDIM x JDIM x kx) ! PHALF Pressure at half model levels ! (dimensioned IDIM x JDIM x kx+1) ! PSFC Surface pressure field ! (dimensioned IDIM x JDIM) ! TIME time of year (time_type) ! LAT latitudes in radians, dimensioned by (IDIMxJDIM) ! KBOT - OPTIONAL; lowest model level index array ! (dimensioned IDIM x JDIM) !=================================================================== ! Arguments (intent out) integer, intent(out), dimension(:,:,:) :: cldtop,cldbas,icld integer, intent(out), dimension(:,:) :: nclds integer, intent(out), dimension(:,:) :: lhight,lhighb, lmidt, lmidb, llowt real, intent(out), dimension(:,:,:) :: cldamt real, intent(out), dimension(:,:,:) :: tempcld,delp_true integer, intent(in), OPTIONAL, dimension(:,:) :: kbot ! OUTPUT ! ------ ! ICLD marker array of cloud types (at cloud levels) ! (dimensioned IDIM x JDIM x kx) ! CLDTOP index of cloud tops (at cloud levels) ! (dimensioned IDIM x JDIM x kx) ! CLDBAS index of cloud bottoms (at cloud levels) ! (dimensioned IDIM x JDIM x kx) ! CLDAMT cloud amount (fraction) (at cloud levels) ! (dimensioned IDIM x JDIM x kx) ! TEMPCLD cloud layer mean temperature (degrees Kelvin) ! (at cloud levels) ! DELP_TRUE true cloud pressure thickness of distinct cloud layers ! (at cloud levels) ! R_UV fractional amount of ultraviolet radiation ! reflected by the clouds ! R_NIRfractional amount of near inrared radiation ! reflected by the clouds ! AB_UVfractional amount of ultraviolet radiation ! absorbed by the clouds ! AB_NIRfractional amount of near inrared radiation ! absorbed by the clouds ! EM_LWemissivity for the clouds ! NCLDS number of (random overlapping) clouds in column and also ! the current # for clouds to be operating on ! LHIGHT vertical level index upper limit for high cloud tops ! a function of lon and latitude (dimensioned IDIMxJDIM) ! LHIGHB vertical level index lower limit for high cloud bases ! a function of lon and latitude (dimensioned IDIMxJDIM) ! LMIDT vertical level index upper limit for mid cloud tops ! a function of lon and latitude (dimensioned IDIMxJDIM) ! LMIDB vertical level index lower limit for mid cloud bases ! a function of lon and latitude (dimensioned IDIMxJDIM) ! LLOWT vertical level index upper limit for low cloud tops ! a function of lon and latitude (dimensioned IDIMxJDIM) ! ! note: vertical level below which no low cloud bases can exist is ! calculated in routines where needed using namelist inputs ! low_lev_cloud_index and nofog !=================================================================== !======================================================================= ! (Intent local) real , dimension(size(rhum,1),size(rhum,2),size(rhum,3)) :: theta real , dimension(size(rhum,1),size(rhum,2),size(rhum,3)) :: rhumcnv, & camtcnv,camtrh,camtw,camtsh,camtsc,camt real, dimension (size(rhum,1),size(rhum,2)) :: camtsh_mx real, dimension (size(phalf,1),size(phalf,2),size(phalf,3)) :: pnorm integer , dimension(size(rhum,1),size(rhum,2),size(rhum,3)) :: icld_k integer , dimension(size(rhum,1),size(rhum,2),3) :: kthclm, kbhclm integer, dimension (size(rhum,1),size(rhum,2)) :: kmaxshl,kminshl ! horizontal dimensions integer idim,jdim ! bottom level vertical index integer kx ! loop index integer k !----------------------------------------------------------------------- ! THETA Potential temperature ! at full model levels (Deg K) ) ! (dimensioned IDIM x JDIM x kx) ! THETA_E Equivalent potential temperature ! at full model levels (Deg K) ) ! (dimensioned IDIM x JDIM x kx) ! RHUMCNV Relative humidity fraction modified for convective clouds ! (dimensioned IDIM x JDIM x kx) ! CAMTCNV tentative deep convective cloud amounts ! (dimensioned IDIM x JDIM x kx) ! CAMTRH tentative rel. humidity cloud amounts ! (dimensioned IDIM x JDIM x kx) ! CAMTW tentative cloud amounts omega corrected clouds ! (dimensioned IDIM x JDIM x kx) ! CAMTSH tentative shallow convective cloud amounts ! (dimensioned IDIM x JDIM x kx) ! CAMTSC tentativemarine stratus cloud amounts ! (dimensioned IDIM x JDIM x kx) ! CAMTSL tentative Slingo marine stratus cloud amounts ! (dimensioned IDIM x JDIM x kx) ! CAMTTL tentative Tim Li marine stratus cloud amounts ! (dimensioned IDIM x JDIM x kx) ! CAMTSH_MX maximum value of shallow convective cloud amount within ! shallow convective layer. (dimensioned IDIM x JDIM) ! CAMT tentative merged cloud amounts ! (dimensioned IDIM x JDIM x kx) ! PNORM Normalized pressure at half or full model levels ! (dimensioned IDIM x JDIM x kx+1) ! ICLD_K tentative merged cloud marker array ! (dimensioned IDIM x JDIM x kx) ! KTHCLM,KBHCLM gcm vert coord indices of climo high, mid, low cloud ! tops and bases, (dimensioned IDIM x JDIM x 3) ! KMINSHL,KMAXSHL indices of vertical levels corresponding to ! final top and base of shallow convective cloud layer ! (dimensioned IDIM x JDIM) !----------------------------------------------------------------------- ! Calculate gcm vetical coord. indices for high, mid and low clouds based ! on observed mean climatology. They will be a function of latitude and ! time of year. Optionally, climatological zonal mean cloud amounts may ! also be returned, presumably for a special "cold" start situation. ! ***** caveats: the current climo cloud setup uses climo average pressure ! ***** values for locating single layer high and middle clouds ! ***** and the bottom and top of climo average low clouds ! ***** in the future a better approach to locating a vertical ! ***** layer range for each type of cloud will be needed ! ***** camt from getcld is only used to complete the argument list ! need to normalize pressures since table is based on assumption that ! surface pressure = 101325 Pa idim = size(rhum,1) jdim = size(rhum,2) kx = size(rhum,3) ! for clouds at half levels ! do k=1,kx+1 ! pnorm(:,:,k)=101325.*phalf(:,:,k)/phalf(:,:,kx+1) ! enddo ! call getcld (time, lat, phalf, kthclm, kbhclm, camt) ! for clouds at full levels do k=1,kx pnorm(:,:,k)=101325.*pfull(:,:,k)/phalf(:,:,kx+1) enddo pnorm(:,:,kx+1) = 101325. ! call getcld (time, lat, pfull, kthclm, kbhclm, camt) call getcld (time, lat, pnorm, kthclm, kbhclm, camt) ! re-initialize camt, icld_k, nclds, cldbas,cldtop camt(:,:,:) = 0.0 icld_k(:,:,:) = 0 cldbas(:,:,:) = 0 cldtop(:,:,:) = 0 nclds(:,:) = 0 lhight(:,:) = high_lev_cloud_index lhighb(:,:) = kthclm(:,:,2) - 1 lmidt(:,:) = kthclm(:,:,2) ! if climo mid cloud top level = climo high cloud base level, reset high base where (kthclm(:,:,2) .eq. kbhclm(:,:,1)) lhighb(:,:) = kbhclm(:,:,1) lmidt(:,:) = kbhclm(:,:,1) + 1 endwhere lmidb(:,:) = kthclm(:,:,3) - 1 llowt(:,:) = kthclm(:,:,3) ! if climo low cloud top level = climo mid cloud base level, reset mid base where (kthclm(:,:,3) .eq. kbhclm(:,:,2)) lmidb(:,:) = kbhclm(:,:,2) llowt(:,:) = kbhclm(:,:,2) + 1 endwhere ! print *, ' climatological cloud tops & bases' ! print *, ' kthclm = ', kthclm ! print *, ' kbhclm = ', kbhclm !----------------------------------------------------------------------- ! calculate potential temperature for use in computing marine stratus and/or ! shallow convective clouds if (linvers .or. l_theqv) then theta(:,:,:) = temp(:,:,:)*(p00/pfull(:,:,:)**Kappa) endif ! calculate deep convective clouds if (lcnvcld) then call cloud_cnv (rhum,cnvcntq,convprc, pfull, phalf, camtcnv, rhumcnv ) endif ! calculate rel humidity clouds if (lrhcnv .and. lcnvcld) then call cloud_rhum (rhumcnv, pnorm, camtrh) else call cloud_rhum (rhum, pnorm, camtrh) endif ! calculate omega corrected rel humidity clouds if (lomega) then call cloud_omga (camtrh, omega, llowt, camtw) endif ! calculate marine stratus clouds ! ******* not implemented at this time ******* if (linvers) then !!! call cloud_m_stratus (...,camtsc ) endif kx = size(rhum,3) ! calculate shallow convective clouds if (l_theqv) then ! need to call shallow_conv_init to set up constants for mylcl call cloud_shallow_conv (theta,omega,pfull,phalf,temp,qmix,camtrh, & camtsh,camtsh_mx,kminshl,kmaxshl,kbot) endif ! merge stratiform cloud types call merge_strat_clouds (camtrh,camtw,camtsc,llowt,camt,icld_k) ! group clouds into high, middle and low call def_hi_mid_low (phalf,lhight,lhighb,lmidt,lmidb,llowt,camt,icld_k) ! merge convective cloud types if (lcnvcld .or. l_theqv) then call merge_cnv_clouds (camtcnv,camtsh,camtsh_mx,kminshl,kmaxshl, & lhight,lhighb,lmidt,lmidb,llowt,camt,icld_k) endif ! calculate vertical indices and tot number of distinct cloud layers ! as a function of lat and longitude. call layr_top_base (icld_k,nclds,cldtop,cldbas,icld) ! calculate vertical, mass weighted cloud amount fraction of distinct cloud ! layers call cldamt_mn (nclds,cldtop,cldbas,phalf,camt,cldamt) ! compute total cloud amount from cloud amount of in distinct cloud layers ! ***** this diagnostice routine is postpooned for now ******* ! call total_cld_amt() ! Define the occurence of anvil cirrus clouds call anvil_cirrus (lhighb,lmidb,nclds,cldtop,cldbas,cnvcntq,lgscldelq, & pfull,phalf,icld) ! Compute the mass weighted mean cloud temperature and pressure thickness ! of each distinct cloud layer call cld_layr_mn_temp_delp (nclds,cldtop,cldbas,temp,phalf,tempcld,delp_true) ! Compute the cum mean water vapor mixing ratio of each distinct cloud layer ! (This is proposed for the future) ! call cld_lay_cum_qmix() !----------------------------------------------------------------------- end subroutine CLDTIM !======================================================================= !############################################################################# subroutine CLOUD_CNV (rhum,cnvcntq,convprc, pfull, phalf, camtcnv, rhumcnv ) !------------------------------------------------------------------- ! This subroutine calculates deep convective cloud amounts !------------------------------------------------------------------- !=================================================================== ! Arguments (intent in) !------------------------------------------------------------------- ! Namelist variables used in this routine (defined at top of module) ! HIGH_LEV_CLOUD_INDEX - level above which no clouds are allowed to form ! (model level index) ! NOFOG - logical switch for not allowing rhum clouds (or fog) ! to occur beneath a certain level (low_lev_cloud_index) -> ! nofog = true ! to allow clouds at the lowest model level -> nofog = false ! LOW_LEV_CLOUD_INDEX - level below which no clouds are allowed to occur ! when nofog = true (model level index) !------------------------------------------------------------------- !------------------------------------------------------------------- ! parameters used in cloud_cnv real, parameter :: beta = .125, cprc_thres = 0.14 real, parameter :: cmax = 0.8, ctower = 0.25 real, parameter :: cnvdelqmn = 0.0 ! BETA - convective cloud regression coefficients ! CPRC_THRES - minimum convective precip amount imposed for ! for work array. This allows the use of the log ! function in the convective cloud regression relation. ! CMAX, CTOWER - factors limiting convective cloud amount ! and vertical profile ! CNVDELQMN - min threshold value of accumulated count of mixing ratio ! change due to convection above which convective clouds are ! allowed ! ( It is set = 0.0 as a parameter, because it not anticipated ! that it would be changed.) !------------------------------------------------------------------- real, intent(in), dimension (:,:,:) :: rhum, cnvcntq, pfull, phalf real, intent(in), dimension (:,:) :: convprc !------------------------------------------------------------------- ! RHUM Relative humidity fraction at full model levels ! (dimensioned IDIM x JDIM x kx) ! CNVCNTQ Accumulated count of change in mix ratio due to conv precip ! at full model levels ! (dimensioned IDIM x JDIM x kx) ! convprc Accumulated conv precip rate summed over all ! full model levels (mm/day ) ! (dimensioned IDIM x JDIM) ! PFULL Pressure at full model levels ! (dimensioned IDIM x JDIM x kx) ! PHALF Pressure at half model levels ! (dimensioned IDIM x JDIM x kx+1) !=================================================================== ! Arguments (intent out) real, intent(out), dimension (:,:,:) :: camtcnv, rhumcnv ! CAMTCNV - convective cloud amount at full model levels ! (dimensioned IDIM x JDIM x kx) ! RHUMCNV - relative humidity adjusted for convective clouds !======================================================================= ! (Intent local) real, dimension (size(rhum,1),size(rhum,2),size(rhum,3)) :: delp real, dimension (size(rhum,1),size(rhum,2)) :: wrkcprc ! DELP pressure thickness of model layers ! (dimensioned IDIM x JDIM x kx) ! DCPRC - weighted convective precip as a function of vertical layer ! (dimensioned IDIM x JDIM x kx) ! WRKCPRC - convective precip work array ! (dimensioned IDIM x JDIM) integer i, j, idim, jdim, kx, lk, lkxm1, kcbtop real alpha !----------------------------------------------------------------------- ! type loop index variable integer k !=================================================================== ! Initialize convective cloud amount array camtcnv = 0.0 ! Define regression coefficient alpha base on coeficient beta, defined ! above as a parameter alpha = -beta*log(cprc_thres) !----------------------------------------------------------------------- ! <><><><><><><><> set up vertical index range <><><><><><><><> !----------------------------------------------------------------------- idim = SIZE(rhum,1) jdim = SIZE(rhum,2) kx = size(rhum,3) ! define lower limit of cloud bottoms if (nofog) then lk = low_lev_cloud_index else lk = kx endif ! no convective clouds allowed in lowest layer lkxm1 = min(lk,kx-1) ! define upper limit of cloud tops kcbtop = high_lev_cloud_index ! calculate pressure thickness of model layers, for vertcal averaging do k=1,kx delp (:,:,k) = phalf(:,:,k+1)-phalf(:,:,k) end do !----------------------------------------------------------------------- ! <><><><><><><><> code to identify convecting levels <><><><><><><><> !----------------------------------------------------------------------- ! use total convective precip amount passed in as argument wrkcprc(:,:) = convprc(:,:) !----------------------------------------------------------------------- ! <><><><><><><><> code to compute conv cloud fraction <><><><><><><><> !----------------------------------------------------------------------- do k=kcbtop,lkxm1 ! check that accumulated count of convective mix ratio changes is greater ! than a minimum threshold (usually 0 ) where (cnvcntq (:,:,k) > cnvdelqmn .and. wrkcprc(:,:) > cprc_thres) camtcnv(:,:,k) = alpha + beta * log(wrkcprc(:,:) ) end where end do ! Impose contraints on convective cloud camtcnv(:,:,:) = min(camtcnv(:,:,:),cmax) do k=kcbtop,lkxm1 do j=1,jdim do i=1,idim if (camtcnv(i,j,k) .lt. 0.0) then print *, ' pe,i,j,k,camtcnv = ', mpp_pe(),i,j,k, & camtcnv(i,j,k) call error_mesg ('cloud_cnv','cloud amount < 0' ,FATAL) endif end do end do end do where (camtcnv(:,:,:) > 0.0 .and. pfull(:,:,:) < pshallow) camtcnv(:,:,:) = ctower * camtcnv(:,:,:) endwhere ! calculate relative humidity adjusted for convective clouds rhumcnv (:,:,:) = (1.0 - camtcnv(:,:,:)) * rhum(:,:,:) !----------------------------------------------------------------------- end subroutine CLOUD_CNV !############################################################################# subroutine CLOUD_RHUM (rhum,pnorm, camtrh) !------------------------------------------------------------------- ! calculates stratiform cloud amounts based on relative humidities !------------------------------------------------------------------- !=================================================================== ! Arguments (intent in) ! Namelist variables used in this routine (defined at top of module) ! HIGH_LEV_CLOUD_INDEX - level above which no clouds are allowed to form ! (model level index) ! NOFOG - logical switch for not allowing rhum clouds (or fog) ! to occur beneath a certain level (low_lev_cloud_index) -> ! nofog = true ! to allow clouds at the lowest model level -> nofog = false ! LOW_LEV_CLOUD_INDEX - level below which no clouds are allowed to occur ! when nofog = true (model level index) !------------------------------------------------------------------- real, intent(in), dimension (:,:,:) :: rhum,pnorm !------------------------------------------------------------------- ! RHUM Relative humidity fraction at full model levels ! (dimensioned IDIM x JDIM x kx) ! PNORM Normalized pressure at half or full model levels ! (dimensioned IDIM x JDIM x kx+1) !=================================================================== ! Arguments (intent out) real, intent(out), dimension (:,:,:) :: camtrh ! CAMTRH - Rel Humidity cloud amount at full model levels ! (dimensioned IDIM x JDIM x kx) !======================================================================= ! (Intent local) real, dimension (size(rhum,1),size(rhum,2)) :: rhc_work ! real, dimension (size(rhum,1),size(rhum,2),size(rhum,3)) :: rhum2 integer kx, lk, npower !----------------------------------------------------------------------- ! type loop index variable integer k,kk !=================================================================== kx = size(rhum,3) ! define cloud amt - rel hum relation as linear or quadratic if (lquadra) then npower = 2 else npower = 1 endif ! if (mpp_pe() == 0) then ! print *, 'npower', npower ! endif ! define lower limit rel hum clouds if (nofog) then lk = low_lev_cloud_index else lk = kx endif ! Initialize rhum cloud amount array camtrh = 0.0 ! !----------------------------------------------------------------------- ! <><><><><><><><> code to compute rhum cloud fraction <><><><><><><><> !----------------------------------------------------------------------- do k=high_lev_cloud_index,lk where (pnorm(:,:,k) .ge. pbounds(nrhc-1)) rhc_work(:,:) = rhc(nrhc) end where do kk=nrhc-1,2 where ((pnorm(:,:,k).lt.pbounds(kk)) .and. & (pnorm(:,:,k).ge.pbounds(kk-1))) rhc_work(:,:) = rhc(kk) end where end do where (pnorm(:,:,k).lt.pbounds(1)) rhc_work(:,:) = rhc(1) end where !BUGFIX ?? ! where (rhum(:,:,k) == rhc_work(:,:)) ! rhum2(:,:,k) = rhum(:,:,k) - 0.00007 !elsewhere ! rhum2(:,:,k) = rhum(:,:,k) ! endwhere where (rhum (:,:,k) > rhc_work(:,:)) camtrh(:,:,k) = min(1.0, ( & (rhum (:,:,k) - rhc_work(:,:))/(1-rhc_work(:,:)) )** npower) elsewhere camtrh(:,:,k) = 0.0 end where end do end subroutine CLOUD_RHUM !############################################################################# subroutine CLOUD_OMGA (camtrh, omega, llowt, camtw) !------------------------------------------------------------------- ! calculates omega corrected cloud amounts ! This subroutine reduces rel hum clouds in regions of decending motion !------------------------------------------------------------------- !=================================================================== ! Arguments (intent in) ! Namelist variables used in this routine (defined at top of module) ! NOFOG - logical switch for not allowing rhum clouds (or fog) ! to occur beneath a certain level (low_lev_cloud_index) -> ! nofog = true ! to allow clouds at the lowest model level -> nofog = false ! LOW_LEV_CLOUD_INDEX - level below which no clouds are allowed to occur ! when nofog = true (model level index) ! WCUT0 - omega cutoff value for omega cloud depletion factor = 0 ! WCUT1 - omega cutoff value for omega cloud depletion factor = 1 !----------------------------------------------------------------------- real, intent(in), dimension (:,:,:) :: camtrh, omega integer, intent(in), dimension (:,:) :: llowt !----------------------------------------------------------------------- ! CAMTRH tentative rel. humidity cloud amounts ! (dimensioned IDIM x JDIM x kx) ! OMEGA Pressure vertical velocity at full model levels ! (dimensioned IDIM x JDIM x kx) ! LLOWT vertical level index upper limit for low cloud tops ! a function of lon and latitude (dimensioned IDIMxJDIM) !----------------------------------------------------------------------- !=================================================================== ! Arguments (intent out) real, intent(out), dimension (:,:,:) :: camtw ! CAMTW - omega cloud amount at full model levels ! (dimensioned IDIM x JDIM x kx) !======================================================================= ! (Intent local) !----------------------------------------------------------------------- ! type loop limits integer kx, lk, idim, jdim integer low_top !----------------------------------------------------------------------- ! type loop index variables integer i,j, k !=================================================================== kx = size(camtrh,3) jdim = size(camtrh,2) idim = size(camtrh,1) ! Initialize omega cloud amount array camtw = 0.0 ! define lower limit for calculating omega corrected clouds if (nofog) then lk = low_lev_cloud_index else lk = kx endif ! llowt is used to define an upper bound for calculating omega corrected ! clouds, i.e., estimate of highest level of a low cloud top !----------------------------------------------------------------------- !=================================================================== !----------------------------------------------------------------------- ! <><><><><><><> code to compute omega corrected cloud fraction <><><> !----------------------------------------------------------------------- do j=1,jdim do i=1,idim low_top = llowt(i,j) do k=low_top,lk ! For regions of upward motion (or zero) set = rhum clouds if (omega (i,j,k) .le. wcut1 ) then camtw(i,j,k) = camtrh(i,j,k) endif ! For regions of descending motion reduce clouds linearly to 0 where ! descent = 3.6mb/hr = .1 n m-2 s-1 if (omega (i,j,k) .gt. wcut1 .and. omega (i,j,k) .lt. wcut0 ) then camtw(i,j,k) = camtrh(i,j,k)*(omega(i,j,k) - wcut0)/(wcut1-wcut0) endif ! For regions of descent >= 3.6mb/hr set to 0 if (omega (i,j,k) .ge. wcut0 ) then camtw(i,j,k) = 0.0 endif end do end do end do !----------------------------------------------------------------------- end subroutine CLOUD_OMGA !############################################################################# subroutine CLOUD_SHALLOW_CONV (theta,omega,pfull,phalf,temp,qmix,camtrh, & camtsh,camtsh_mx,kminshl,kmaxshl,kbot) !------------------------------------------------------------------- ! This subroutine calculates shallow convective cloud amounts !------------------------------------------------------------------- !=================================================================== ! Arguments (intent in) ! Namelist variables used in this routine (defined at top of module) ! PSHALLOW - top of shallow convective layer (pressure level - n/m**2 ) ! L_THEQV - logical switch for turning on calculation of shallow convective ! clouds - true for shallow convective clouds, otherwise false ! WCUT1 - omega cutoff value for omega cloud depletion factor = 1 ! LOMEGA - logical switch for turning on omega correction to rhum ! clouds - true for omega correction, otherwise false ! LINVERS - logical switch for turning on calculation of marine stratus ! clouds - true for marine stratus, otherwise false ! HIGH_LEV_CLOUD_INDEX - level above which no clouds are allowed to form ! (model level index) !------------------------------------------------------------------- ! parameters used in cloud_shallow_conv real, parameter :: twcfac = 0.2, crtcons = 0.0 ! TWCFAC - scaling factor for computing shallow conv cloud amt (= 0.2) ! CRTCONS - crit value of d(Theta_e)/dP (= 0.0) real, intent(in), dimension (:,:,:) :: & theta,omega,pfull,phalf,temp,qmix,camtrh integer, intent(in), OPTIONAL, dimension(:,:) :: kbot ! THETA Potential temperature at full model levels ( Deg K ) ! (dimensioned IDIM x JDIM x kx) ! OMEGA Pressure vertical velocity at full model levels ! (dimensioned IDIM x JDIM x kx) ! PFULL Pressure at full model levels ! (dimensioned IDIM x JDIM x kx) ! PHALF Pressure at half model levels ! (dimensioned IDIM x JDIM x kx+1) ! TEMP Temperature (Deg K) at full model levels ! (dimensioned IDIM x JDIM x kx) ! QMIX Mixing Ratio at full model levels ! (dimensioned IDIM x JDIM x kx) ! CAMTRH tentative rel. humidity cloud amounts ! (dimensioned IDIM x JDIM x kx) ! KBOT OPTIONAL; lowest model level index array ! (dimensioned IDIM x JDIM) !=================================================================== ! Arguments (intent out) real, intent(out), dimension (:,:,:) :: camtsh real, intent(out), dimension (:,:) :: camtsh_mx integer, intent(out), dimension (:,:) :: kminshl, kmaxshl ! CAMTSH - shallow convective cloud amount at full model levels ! (dimensioned IDIM x JDIM x kx) ! CAMTSH_MX maximum value of shallow convective cloud amount within ! shallow convective layer. (dimensioned IDIM x JDIM) ! KMINSHL,KMAXSHL indices of vertical levels corresponding to ! final top and base of shallow convective cloud layer ! (dimensioned IDIM x JDIM) !======================================================================= ! (Intent local) real, dimension(size(temp,1),size(temp,2),size(temp,3)) :: & qsat,dcldnrm,theta_e,dtheta_e real, dimension (size(temp,1),size(temp,2)) :: plcl integer, dimension (size(temp,1),size(temp,2)) :: ksiglcl,kshallow ! QSAT - saturation vapor pressure ! (dimensioned IDIM x JDIM x KX) ! PLCL - pressure of the lifting condensation level (LCL) ! (dimensioned IDIM x JDIM) ! KSIGLCL - index of nearest vertical level to the LCL ! (dimensioned IDIM x JDIM) ! DCLDNRM - normalization factor ! (dimensioned IDIM x JDIM x KX) ! THETA_E Equivalent potential temperature at full model levels ( Deg K ) ! (dimensioned IDIM x JDIM x kx) ! DTHETA_E - vertical (p) derivative of theta_e ! (dimensioned IDIM x JDIM x KX) ! define 2-D work arrays real, dimension(size(temp,1),size(temp,2)) :: xy1, xy2, xy3 !----------------------------------------------------------------------- ! type local variables integer i, j, k, kmin, kmax, kmin0, kmax0, idim, jdim,lk, lkxm1, kx, & kcbtop, kcount, knt_shl !======================================================================= !======================================================================= idim = SIZE(Temp,1) jdim = SIZE(Temp,2) kx = SIZE(Temp,3) ! define lower limit of cloud bottom if (nofog) then lk = low_lev_cloud_index else lk = kx endif ! no convective clouds allowed in lowest layer lkxm1 = min(lk,kx-1) ! define upper limit of cloud tops kcbtop = high_lev_cloud_index ! find index of vertical level closest to but not above the ! shallow convective lid (pshallow) kshallow = 99999 do k = kx,1,-1 where (pfull(:,:,k) .ge. pshallow) kshallow(:,:) = k end where end do ! Initialize shallow convective cloud amount and other related arrays camtsh = 0.0 camtsh_mx = 0.0 dcldnrm = 0.0 ! calculate lcl ! --- saturation mixing ratio call compute_qs (Temp, pfull, qsat) ! calculate equivalent potential temperature theta_e(:,:,:) = theta(:,:,:)*exp(HLv*qmix(:,:,:)/ & (Cp_Air*temp(:,:,:))) !======================================================================= ! --- CALCULATE THE LIFTING CONDENSATION LEVEL, IE CLOUB BASE !======================================================================= ! The optional argument array kbot allows for use of discontinuous vertical ! coordinate systems. However, if the kbot index value relates to a level ! which is well above the local surface, the resulting lcl may not be ! physically correct. if( PRESENT( kbot ) ) then do j=1,jdim do i=1,idim k = kbot(i,j) xy1(i,j) = Temp(i,j,k) xy2(i,j) = MIN(MAX( qmix(i,j,k), 1.0e-6), qsat(i,j,k) ) xy3(i,j) = pfull(i,j,k) end do end do else xy1(:,:) = Temp(:,:,kx) xy2(:,:) = MIN(MAX( qmix(:,:,kx), 1.0e-6), qsat(:,:,kx) ) xy3(:,:) = pfull(:,:,kx) end if CALL MYLCL( xy1, xy2, xy3, phalf, plcl, ksiglcl ) !======================================================================= ! --- Calculate d(theta_e)/dp !======================================================================= ! set dtheta_e = 0 at the first level dtheta_e(:,:,1) = 0.0 do k=2,kx-1 dtheta_e(:,:,k) = & 0.5*(theta_e(:,:,k+1)-theta_e(:,:,k))/(pfull(:,:,k+1)-pfull(:,:,k)) + & 0.5*(theta_e(:,:,k)-theta_e(:,:,k-1))/(pfull(:,:,k)-pfull(:,:,k-1)) end do ! set dtheta_e at bottom level = detheta_e at next level up ! (note: bottom level value is not used in current scheme) dtheta_e(:,:,kx) = dtheta_e(:,:,kx-1) !======================================================================= ! --- Shallow Convective Cloud Calculation !======================================================================= kmaxshl(:,:) = min(lkxm1,ksiglcl(:,:)) ! tentatively set kminshl to kshallow kminshl(:,:) = kshallow(:,:) do k = lkxm1,kcbtop,-1 where (((k.ge.kshallow(:,:)) .and. (k.le.kmaxshl(:,:))) & .and. (dtheta_e(:,:,k).ge.crtcons ) & .and. (pfull(:,:,k) .le. plcl(:,:) ) ) dcldnrm(:,:,k) = twcfac end where end do ! Apply constraint that the shallow convective layer is confined to the ! lowest contigous layer between plcl and pshallow do j=1,jdim do i=1,idim ! no shallow conv clouds allowed in regions where lowest full P level ! is above the specified shallow convective lid if (kshallow(i,j) .le. kx) then kmax=kmaxshl(i,j) kmin=kshallow(i,j) kcount = 0 knt_shl = 0 do k = kmax,kmin,-1 if (pfull(i,j,k) .le. plcl(i,j) ) then if (dcldnrm(i,j,k).eq.0.0 .and. knt_shl.ge.1) then dcldnrm(i,j,k-1) = 0.0 kcount = kcount + 1 ! kminshl calculation is completed here if (kcount .eq. 1) then kminshl(i,j) = k endif else if (dcldnrm(i,j,k) .gt. 0.0) knt_shl = knt_shl + 1 endif endif end do else dcldnrm(i,j,:) = 0.0 endif end do end do ! calculate cloud amounts do k = lkxm1,kcbtop,-1 where ((k.ge.kminshl(:,:)) .and. (k .le.kmaxshl(:,:)) & .and. dcldnrm(:,:,k).gt.0.0 .and. (omega(:,:,k).gt.wcut1) ) camtsh(:,:,k) = dcldnrm(:,:,k)*camtrh(:,:,k) camtsh_mx (:,:) = max(camtsh_mx (:,:),camtsh (:,:,k) ) end where end do ! impose maximum cloud overlap within the contiguous shallow convective ! cloud layer. Also check that an rhum cloud exists and that the condition ! on vertical motion is satisfied do k = lkxm1,kcbtop,-1 where ((k.ge.kminshl(:,:)) .and. (k .le.kmaxshl(:,:)) & .and. dcldnrm(:,:,k).gt.0.0 .and. (camtrh(:,:,k).gt.0.0) & .and. (omega(:,:,k).gt.wcut1) ) camtsh(:,:,k) = camtsh_mx (:,:) elsewhere camtsh(:,:,k) = 0.0 endwhere end do ! redefine kminshl and kmaxshl based on actual shallow conv clouds do j=1,jdim do i=1,idim if (camtsh_mx(i,j) .gt. 0.0) then kmax = kmaxshl(i,j) kmax0 = kmaxshl(i,j) kmin = kminshl(i,j) kmin0 = kminshl(i,j) ! find actual shallow conv cloud base by searching upwards do k=kmax0,kmin0,-1 if (camtsh(i,j,k) .gt. 0.0) then kmax = k go to 100 endif end do 100 continue ! find actual shallow conv cloud top by searching downward do k=kmin0,kmax if (camtsh(i,j,k) .gt. 0.0) then kmin = k go to 200 endif end do 200 continue kmaxshl(i,j) = kmax kminshl(i,j) = kmin endif end do end do !----------------------------------------------------------------------- !----------------------------------------------------------------------- end subroutine CLOUD_SHALLOW_CONV !####################################################################### subroutine MERGE_STRAT_CLOUDS (camtrh,camtw,camtsc,llowt,camt,icld_k) !------------------------------------------------------------------- ! This subroutine determines a dominant stratiform cloud amount !------------------------------------------------------------------- !=================================================================== ! Arguments (intent in) ! Namelist variables used in this routine (defined at top of module) ! LOMEGA - logical switch for turning on omega correction to rhum ! clouds - true for omega correction, otherwise false ! LINVERS - logical switch for turning on calculation of marine stratus ! clouds - true for marine stratus, otherwise false ! NOFOG - logical switch for not allowing rhum clouds (or fog) ! to occur beneath a certain level (low_lev_cloud_index) -> ! nofog = true ! to allow clouds at the lowest model level -> nofog = false ! LOW_LEV_CLOUD_INDEX - level below which no clouds are allowed to occur ! when nofog = true (model level index) real, intent(in), dimension (:,:,:) :: camtrh, camtw, camtsc integer, intent(in), dimension (:,:) :: llowt ! CAMTRH tentative cloud amounts from stratiform clouds ! (dimensioned IDIM x JDIM x kx) ! CAMTW tentative cloud amounts omega corrected clouds ! (dimensioned IDIM x JDIM x kx) ! CAMTSC tentativemarine stratus cloud amounts ! (dimensioned IDIM x JDIM x kx) ! LLOWT vertical level index upper limit for low cloud tops ! a function of lon and latitude (dimensioned IDIMxJDIM) !=================================================================== ! Arguments (intent out) real, intent(out), dimension (:,:,:) :: camt integer, intent(out), dimension (:,:,:) :: icld_k ! CAMT tentative (only stratiform) merged cloud amounts ! (dimensioned IDIM x JDIM x kx) ! ICLD_K tentative merged cloud marker array ! (dimensioned IDIM x JDIM x kx) !======================================================================= ! (Intent local) integer kx, lk, idim, jdim, i, j, k integer low_top !======================================================================= !======================================================================= idim = SIZE(camtrh,1) jdim = SIZE(camtrh,2) kx = SIZE(camtrh,3) ! define lower limit for stratiform clouds if (nofog) then lk = low_lev_cloud_index else lk = kx endif ! Initialize output cloud amounts to rel hum cloud amounts and define ! cloud type marker values camt(:,:,:) = camtrh(:,:,:) where ( camtrh(:,:,:) .gt. 0.0) icld_k(:,:,:) = 1 endwhere ! assign omega corrected cloud amounts if (lomega) then do j=1,jdim do i=1,idim low_top =llowt(i,j) do k=low_top,lk if ( camtw(i,j,k).gt. 0.0) then camt(i,j,k) = camtw(i,j,k) if ( camtw(i,j,k) .lt. camtrh(i,j,k) ) then icld_k(i,j,k) = 2 endif endif end do end do end do endif ! assign marine stratus cloud amounts if (linvers) then where ( camtsc(:,:,:) .gt. camt(:,:,:) ) camt(:,:,:) = camtsc(:,:,:) icld_k(:,:,:) = 3 endwhere endif !----------------------------------------------------------------------- end subroutine MERGE_STRAT_CLOUDS !####################################################################### subroutine DEF_HI_MID_LOW (phalf,lhight,lhighb,lmidt,lmidb,llowt,camt,icld_k) !------------------------------------------------------------------- ! This subroutine transforms a vertical profile of stratiform cloud ! amounts into "high", "middle" and "low" cloud layers ! if a "thick" cloud option is in effect the current implementation ! extends the cloud up one layer (with current tentative cloud amounts) !------------------------------------------------------------------- !=================================================================== ! Arguments (intent in) ! Namelist variables used in this routine (defined at top of module) ! LTHICK_HIGH - logical variable = true -> allow possibility of raising ! high cloud tops one sigma level to increase their thickness ! from 1 to 2 sigma levels; otherwise they remain thin ! (1 sigma level) ! LTHICK_MID - logical variable = true -> allow possibility of raising ! mid cloud tops one sigma level to increase their thickness ! from 1 to 2 sigma levels; otherwise they remain thin ! (1 sigma level) ! LTHICK_LOW - logical variable = true -> allow possibility of raising ! low cloud tops one sigma level to increase their thickness ! from 1 to 2 sigma levels; otherwise they remain thin ! (1 sigma level) ! NOFOG - logical switch for not allowing rhum clouds (or fog) ! to occur beneath a certain level (low_lev_cloud_index) -> ! nofog = true ! to allow clouds at the lowest model level -> nofog = false ! LOW_LEV_CLOUD_INDEX - level below which no clouds are allowed to occur ! when nofog = true (model level index) real, intent(in), dimension (:,:,:) :: phalf integer,intent(in), dimension(:,:) :: lhight, lhighb, lmidt, lmidb, llowt ! PHALF Pressure at half model levels ! (dimensioned IDIM x JDIM x kx+1) ! LHIGHT vertical level index upper limit for high cloud tops ! a function of lon and latitude (dimensioned IDIMxJDIM) ! LHIGHB vertical level index lower limit for high cloud bases ! a function of lon and latitude (dimensioned IDIMxJDIM) ! LMIDT vertical level index upper limit for mid cloud tops ! a function of lon and latitude (dimensioned IDIMxJDIM) ! LMIDB vertical level index lower limit for mid cloud bases ! a function of lon and latitude (dimensioned IDIMxJDIM) ! LLOWT vertical level index upper limit for low cloud tops ! a function of lon and latitude (dimensioned IDIMxJDIM) ! ! !=================================================================== ! Arguments (intent inout) real, intent(inout), dimension (:,:,:) :: camt integer, intent(inout), dimension (:,:,:) :: icld_k ! CAMT tentative merged cloud amounts ! (dimensioned IDIM x JDIM x kx) ! ICLD_K tentative merged cloud marker array ! (dimensioned IDIM x JDIM x kx) ! !======================================================================= ! (Intent local) ! parameter used in def_hi_mid_low real, parameter :: trshld_camt = 0.25 integer, parameter :: nmax = 2 ! TRSHLD_CAMT - This is a cloud amount threshold value, used in ! conjunction with the thick cloud namelist options. If a thick ! cloud option is on then the level above the cloud amount max ! for a particular cloud type is considered to be part of an ! extended cloud if it is within this fraction of the max cloud amount. ! real, dimension(SIZE(camt,1),SIZE(camt,2),SIZE(camt,3)) :: delp,c_work integer, dimension(SIZE(camt,1),SIZE(camt,2),SIZE(camt,3)) :: ic_work integer, dimension(1) :: kindex ! DELP pressure thickness of model layers ! (dimensioned IDIM x JDIM x kx) ! C_WORK cloud amount work array ! (dimensioned IDIM x JDIM x kx) ! IC_WORK cloud marker work array ! (dimensioned IDIM x JDIM x kx) ! NMAX maximum number of model levels allowed in a thick cloud ! loop limits and indices integer idim,jdim,kx, lk, i, j, k, kbtm, ktop ! special setting of max levels for special cases of mid & high thick clouds integer nmax_mh ! special threshold value for allowing thick middle and high clouds ! when cloud amounts at contiguous levels are the same value (typically 1.0) real trshld_mh idim = SIZE(camt,1) jdim = SIZE(camt,2) kx = SIZE(camt,3) !======================================================================= trshld_mh = .0001 nmax_mh = kx !======================================================================= ! define lower limit for stratiform clouds if (nofog) then lk = low_lev_cloud_index else lk = kx endif !----------------------------------------------------------------------- ! caclulate stratiform l-m-h cloud amounts for thin or thick clouds ! current implementation of thick clouds will extend max stratiform ! cloud layer upward 1 level - assuming the level ! above the max level is within a threshold cloud amount ! [note: the maxloc function used below will choose the lowest value index ! in case of two identical max values. In the extremely unlikely event of ! two identical max cloud amounts in a particular column section, the ! maxloc selection priority will result in a bias toward higher clouds.] c_work(:,:,:) = 0.0 ic_work(:,:,:) = 0 ! calculate pressure thickness of model layers, for vertcal averaging do k=1,kx delp (:,:,k) = phalf(:,:,k+1)-phalf(:,:,k) end do do j = 1,jdim do i = 1,idim ! low clouds kbtm = min(lk,kx-1) ktop = llowt(i,j) kindex(:) = maxloc (camt(i,j,ktop:kbtm)) k=kindex(1) + ktop - 1 c_work(i,j,k) = camt(i,j,k) ! for low clouds add 300 to previous cloud marker type ic_work(i,j,k) = icld_k(i,j,k) + 300 if (lthick_low) then ! do k = kbtm,ktop,-1 call THICK_CLOUDS (camt,delp,trshld_camt,i,j,k,ktop,kbtm,nmax,c_work,ic_work) ! end do endif ! mid clouds kbtm = lmidb(i,j) ktop = lmidt(i,j) kindex(:) = maxloc (camt(i,j,ktop:kbtm)) k=kindex(1) + ktop - 1 c_work(i,j,k) = camt(i,j,k) ! for mid clouds add 200 to previous cloud marker type ic_work(i,j,k) = icld_k(i,j,k) + 200 if (lthick_mid) then call THICK_CLOUDS (camt,delp,trshld_mh,i,j,k,ktop,kbtm,nmax_mh,c_work,ic_work) ! call THICK_CLOUDS (camt,delp,trshld_camt,i,j,k,ktop,kbtm,nmax,c_work,ic_work) endif ! high clouds kbtm = lhighb(i,j) ktop = lhight(i,j) kindex(:) = maxloc (camt(i,j,ktop:kbtm)) k=kindex(1) + ktop - 1 c_work(i,j,k) = camt(i,j,k) ! for high clouds add 100 to previous cloud marker type ic_work(i,j,k) = icld_k(i,j,k) + 100 if (lthick_high) then call THICK_CLOUDS (camt,delp,trshld_mh,i,j,k,ktop,kbtm,nmax_mh,c_work,ic_work) ! call THICK_CLOUDS (camt,delp,trshld_camt,i,j,k,ktop,kbtm,nmax,c_work,ic_work) endif end do end do ! Store the resulting values of c_work and ic_work back into camt and icld_k where (camt(:,:,:) .eq. 0.0) ic_work(:,:,:) = 0 endwhere ! Thick clouds are currently not allowed to extend to the lowest level, but ! fog is allowed there if nofog = f do k=1,kx-1 camt(:,:,k) = c_work(:,:,k) icld_k(:,:,k) = ic_work(:,:,k) end do ! at lowest level check for fog and assign low cloud type where (icld_k(:,:,kx).gt.0) c_work(:,:,kx) = icld_k(:,:,kx) + 300 endwhere icld_k(:,:,kx) = c_work(:,:,kx) end subroutine DEF_HI_MID_LOW !####################################################################### subroutine THICK_CLOUDS (camt,delp,trshld_camt,i,j,kk,ktop,kbtm,nmax, & c_work,ic_work) !------------------------------------------------------------------- ! This subroutine returns "thick" cloud amounts to routine def_hi_mid_low ! if the criteria for those clouds are met. First the level of the max ! cloud amount in a low, middle or high cloud regime is determined. ! A thick cloud layer consisting of a maximum of "nmax" contiguous levels ! is determined to exist by looking for a cloud amount value one level ! above that is within a specified threshold amount (trshld_camt). If ! a contiguous cloud is found and nmax = 2 the search is over. If no ! contiguous cloud is found or nmax > 2 the search continues one level ! below. If contiguous cloud is found both above and below and nmax > 3 ! the search will continue using the same threshold criteria 2 levels ! above. The search will continue to "zigzag" in this manner as long as ! no maximum cloud thickness criteria has not been violated and as long ! as there is no break in contiguous cloud amounts as described above. ! Eventually a pressure thickness value will overide the nmax levels ! criteria for determining the maximum thickness of the clouds. Another ! factor limiting the thickness of clouds is that the cloud extent must ! remain within its respective high, middle, low cloud regime boundaries ! ( i.e., ktop to kbtm). !------------------------------------------------------------------- ! ! !=================================================================== ! Arguments (intent in) real, intent(in), dimension (:,:,:) :: camt,delp real, intent(in) :: trshld_camt integer, intent(in) :: i,j,kk,ktop,kbtm,nmax ! CAMT tentative merged cloud amounts ! (dimensioned IDIM x JDIM x kx) ! DELP pressure thickness of model layers ! (dimensioned IDIM x JDIM x kx) ! I,J,KK indicies for longitude, latitude and ! vertical level (at level of max cloud amount) ! KTOP,KBTM - vertical index valid range for cloud top, cloud bottom ! NMAX - maximum number of model levels allowed in a thick cloud ! TRSHLD_CAMT - This is a cloud amounht threshold value, used in ! conjunction with the thick cloud namelist options. ! If a thick cloud option is on then level(s) adjacent ! to the level of max cloud amount for a particular ! cloud type is considered to be part of an extended cloud ! if it is within this fraction of the max cloud amount. ! ! !======================================================================= ! Arguments (intent inout) real, intent(inout), dimension(:,:,:) :: c_work integer,intent(inout), dimension(:,:,:) :: ic_work ! C_WORK cloud amount work array ! (dimensioned IDIM x JDIM x kx) ! IC_WORK cloud marker work array ! (dimensioned IDIM x JDIM x kx) ! !======================================================================= ! (Intent local) integer ncount,no_up,no_down,k,c_top,c_btm ! loop index integer ki,pass_count real delp_sum,thick_cld_amt ! NCOUNT - counter for number of model levels in a thick cloud ! NO_UP - if = 1 -> no more upward extension of cloud is possible ! NO_DOWN - if = 1 -> no more downward extension of cloud is possible ! K - vertical level index of thick cloud as it being extended ! C_TOP, C_BTM - Top and bottom of calculated thick cloud ! DELP_SUM - sum of pressure weights for thcik cloud layer ! THICK_CLD_AMT - pressure weighted thick cloud amount !======================================================================= !======================================================================= ncount = 1 no_up = 0 no_down = 0 k = kk c_top = kk c_btm = kk pass_count = 0 ! no clouds to process if (camt(i,j,kk) .eq.0.0) go to 300 100 continue pass_count = pass_count + 1 if (no_up .eq. 0) then ! First check to extend cloud upward k = c_top if (camt(i,j,k) .gt.0.0 .and. k.gt.ktop .and. k.le.kbtm) then if ((camt(i,j,kk)-camt(i,j,k-1)).le.trshld_camt) then c_top = c_top - 1 ic_work(i,j,k-1) = ic_work(i,j,k) ncount = ncount + 1 else ! no upward extension of cloud possible no_up = 1 endif else no_up = 1 endif ! check for reaching limiting number of levels in a thick cloud if (ncount .eq. nmax) go to 200 endif ! Next check to extend cloud downward if (no_down .eq. 0) then k = c_btm if (camt(i,j,k) .gt.0.0 .and. k.ge.ktop .and. k.lt.kbtm) then if ((camt(i,j,kk)-camt(i,j,k+1)).le.trshld_camt) then c_btm = c_btm + 1 ic_work(i,j,k+1) = ic_work(i,j,k) ncount = ncount + 1 else ! no downward extension of cloud possible no_down = 1 endif else no_down = 1 endif ! check for reaching limiting number of levels in a thick cloud if (ncount .eq. nmax) go to 200 endif ! Further tests for completion of thick cloud calculations ! Make sure upward or downward extension of the cloud is still allowed ! and the cloud level index is still within the allowable range if ( ( no_up.eq.1 .or. c_top.eq.ktop) .and. & ( no_down.eq.1 .or. c_btm.eq.kbtm) ) go to 200 ! loop back to do another pass at extending the thick cloud go to 100 200 continue ! calculate press thickness weighted thick clouds delp_sum = 0.0 do ki=c_top,c_btm delp_sum = delp_sum + delp(i,j,ki) end do thick_cld_amt = 0.0 do ki=c_top,c_btm thick_cld_amt = thick_cld_amt + camt(i,j,ki)*delp(i,j,ki) end do thick_cld_amt = thick_cld_amt/delp_sum ! broadcast thick cloud amount throughout thick cloud layer do ki=c_top,c_btm c_work(i,j,ki) = thick_cld_amt end do 300 continue !=================================================================== end subroutine THICK_CLOUDS !####################################################################### subroutine MERGE_CNV_CLOUDS (camtcnv,camtsh,camtsh_mx,kminshl,kmaxshl, & lhight,lhighb,lmidt,lmidb,llowt,camt,icld_k) !------------------------------------------------------------------- ! This subroutine adds convective cloud amounts to the already merged ! stratiform cloud amounts !------------------------------------------------------------------- !=================================================================== ! Arguments (intent in) !------------------------------------------------------------------- ! Namelist variables used in this routine (defined at top of module) ! LCNVCLD - logical switch for turning on calculation of deep convective ! clouds - true for deep convective clouds, otherwise false ! L_THEQV - logical switch for turning on calculation of shallow convective ! clouds - true for shallow convective clouds, otherwise false ! HIGH_LEV_CLOUD_INDEX - level above which no clouds are allowed to form ! (model level index) ! NOFOG - logical switch for not allowing rhum clouds (or fog) ! to occur beneath a certain level (low_lev_cloud_index) -> ! nofog = true ! to allow clouds at the lowest model level -> nofog = false ! LOW_LEV_CLOUD_INDEX - level below which no clouds are allowed to occur ! when nofog = true (model level index) !!------------------------------------------------------------------- real, intent(in), dimension (:,:,:) :: camtcnv, camtsh real, intent(in), dimension (:,:) :: camtsh_mx integer, intent(in), dimension (:,:) :: kminshl, kmaxshl integer, intent(in), dimension(:,:) :: lhight, lhighb, lmidt, lmidb, llowt ! CAMTCNV - deep convective cloud amounts at full model levels ! (dimensioned IDIM x JDIM x kx) ! CAMTSH - shallow convective cloud amount at full model levels ! (dimensioned IDIM x JDIM x kx) ! CAMTSH_MX maximum value of shallow convective cloud amount within ! shallow convective layer. (dimensioned IDIM x JDIM) ! KMINSHL,KMAXSHL indices of vertical levels corresponding to ! final top and base of shallow convective cloud layer ! (dimensioned IDIM x JDIM) ! LHIGHT vertical level index upper limit for high cloud tops ! a function of lon and latitude (dimensioned IDIMxJDIM) ! LHIGHB vertical level index lower limit for high cloud bases ! a function of lon and latitude (dimensioned IDIMxJDIM) ! LMIDT vertical level index upper limit for mid cloud tops ! a function of lon and latitude (dimensioned IDIMxJDIM) ! LMIDB vertical level index lower limit for mid cloud bases ! a function of lon and latitude (dimensioned IDIMxJDIM) ! LLOWT vertical level index upper limit for low cloud tops ! a function of lon and latitude (dimensioned IDIMxJDIM) ! !=================================================================== ! Arguments (intent inout) real, intent(inout), dimension (:,:,:) :: camt integer, intent(inout), dimension (:,:,:) :: icld_k ! CAMT tentative merged cloud amounts ! (dimensioned IDIM x JDIM x kx) ! ICLD_K tentative merged cloud marker array ! (dimensioned IDIM x JDIM x kx) !======================================================================= ! (Intent local) real, dimension(SIZE(camt,1),SIZE(camt,2)) :: cmax_strat real, dimension(SIZE(camt,1),SIZE(camt,2),SIZE(camt,3)) :: c_work integer, dimension(SIZE(camt,1),SIZE(camt,2),SIZE(camt,3)) :: ic_work ! CMAX_STRAT max stratiform cloud amount work array ! (dimensioned IDIM x JDIM) ! C_WORK cloud amount work array ! (dimensioned IDIM x JDIM x kx) ! IC_WORK cloud marker work array ! (dimensioned IDIM x JDIM x kx) integer i,j,k,kx, lk, kcbtop, kbtm , ktop, idim, jdim !======================================================================= !======================================================================= idim = SIZE(camt,1) jdim = SIZE(camt,2) kx = SIZE(camt,3) if (nofog) then lk = low_lev_cloud_index else lk = kx endif ! define upper limit of cloud tops kcbtop = high_lev_cloud_index ! initialize work arrays c_work(:,:,:) = camt(:,:,:) ! decide whether or not to replace stratiform low cloud within the ! contiguous shallow convective layer with shallow convective cloud if (l_theqv) then ! find maximum stratiform low cloud amount within shallow conv layer ! assume shallow convective layer is always below upper bound for low clouds do j = 1,jdim do i = 1,idim kbtm = kmaxshl(i,j) ktop = kminshl(i,j) cmax_strat(i,j) = 0.0 do k = kbtm,ktop,-1 cmax_strat(i,j) = max (cmax_strat(i,j),camt(i,j,k)) end do end do end do ! check to see if shallow convective cloud will be dominant type ! if so reset cloud amounts to shallow amount and reset cloud marker do k = lk,kcbtop,-1 where ( (k.ge.kminshl(:,:)) .and. (k.le.kmaxshl(:,:)) .and. & (camtsh_mx(:,:) .gt. cmax_strat(:,:)) ) c_work(:,:,k) = camtsh (:,:,k) icld_k(:,:,k) = 304 endwhere end do endif ! randomly overlap stratiform and deep convective cloud amounts. ! also check to see if the convective cloud is dominant based on ! a comparison of cloud amounts (pre-random overlap), if so assign ! cloud marker values based on H-M-L and cloud type (convective = 5) ic_work(:,:,:) = mod(icld_k(:,:,:),100) if (lcnvcld) then do k = lk,kcbtop,-1 where ( ic_work(:,:,k) .le. 3) c_work(:,:,k) = camtcnv(:,:,k) + (1.0-camtcnv(:,:,k))*camt(:,:,k) endwhere ! low cloud range where ( (camtcnv(:,:,k).gt.camt(:,:,k)) .and. & (k.ge.llowt(:,:) .and. k.le.lk) .and. (ic_work(:,:,k).le.3) ) icld_k(:,:,k) = 305 endwhere ! mid cloud range where ( (camtcnv(:,:,k).gt.camt(:,:,k)) .and. & (k.ge.lmidt(:,:) .and. k.le.lmidb(:,:)) ) icld_k(:,:,k) = 205 endwhere ! high cloud range where ( (camtcnv(:,:,k).gt.camt(:,:,k)) .and. & (k.ge.lhight(:,:) .and. k.le.lhighb(:,:)) ) icld_k(:,:,k) = 105 endwhere end do endif ! re-define convective cloud amount as the max of deep and shallow convective ! cloud at levels where shallow convective clouds are currently defined. if (lcnvcld .and. l_theqv) then where ( (camtcnv(:,:,:).gt.c_work(:,:,:)) .and. & (ic_work(:,:,:) .eq. 4 ) ) c_work(:,:,:) = camtcnv(:,:,:) icld_k(:,:,:) = 305 endwhere endif ! Store the resulting values of c_work back into camt camt(:,:,:) = c_work(:,:,:) !----------------------------------------------------------------------- end subroutine MERGE_CNV_CLOUDS !####################################################################### subroutine LAYR_TOP_BASE (icld_k,nclds,cldtop,cldbas,icld) !------------------------------------------------------------------- ! This subroutine determines cloud layers, tops and bases !------------------------------------------------------------------- !=================================================================== ! Arguments (intent in) ! Namelist variables used in this routine (defined at top of module) ! LTHICK_HIGH - logical variable = true -> allow possibility of raising ! high cloud tops one sigma level to increase their thickness ! from 1 to 2 sigma levels; otherwise they remain thin ! (1 sigma level) ! LTHICK_MID - logical variable = true -> allow possibility of raising ! mid cloud tops one sigma level to increase their thickness ! from 1 to 2 sigma levels; otherwise they remain thin ! (1 sigma level) ! LTHICK_LOW - logical variable = true -> allow possibility of raising ! low cloud tops one sigma level to increase their thickness ! from 1 to 2 sigma levels; otherwise they remain thin ! (1 sigma level) ! HIGH_LEV_CLOUD_INDEX - level above which no clouds are allowed to form ! (model level index) ! NOFOG - logical switch for not allowing rhum clouds (or fog) ! to occur beneath a certain level (low_lev_cloud_index) -> ! nofog = true ! to allow clouds at the lowest model level -> nofog = false ! LOW_LEV_CLOUD_INDEX - level below which no clouds are allowed to occur ! when nofog = true (model level index) !------------------------------------------------------------------- integer, intent(in), dimension (:,:,:) :: icld_k ! ICLD_K tentative merged cloud marker array ! (dimensioned IDIM x JDIM x kx) !=================================================================== ! Arguments (intent out) integer, intent(out), dimension(:,:,:) :: cldtop,cldbas,icld integer, intent(out), dimension(:,:) :: nclds ! NCLDS number of (random overlapping) clouds in column and also ! the current # for clouds to be operating on ! (dimensioned IDIM x JDIM ) ! CLDTOP index of cloud tops (at cloud levels) ! (dimensioned IDIM x JDIM x kx) ! CLDBAS index of cloud bottoms (at cloud levels) ! (dimensioned IDIM x JDIM x kx) ! ICLD marker array of cloud types (at cloud levels) ! (dimensioned IDIM x JDIM x kx) !======================================================================= ! (Intent local) integer, dimension(SIZE(icld_k,1),SIZE(icld_k,2),SIZE(icld_k,3)) :: & ctop_work,cbas_work,ic_work ! CTOP_WORK cloud top work array (at model levels) ! (dimensioned IDIM x JDIM x kx) ! CBAS_WORK cloud bsse work array (at model levels) ! (dimensioned IDIM x JDIM x kx) ! IC_WORK cloud marker work array ! (dimensioned IDIM x JDIM x kx) integer k,kcbtop,i,j,idim, jdim, kx, lk, kp, kpmax, maxcld !======================================================================= !======================================================================= idim = SIZE(icld_k,1) jdim = SIZE(icld_k,2) kx = SIZE(icld_k,3) ! define lower limit of cloud bottom if (nofog) then lk = low_lev_cloud_index else lk = kx endif ! define upper limit of cloud tops kcbtop = high_lev_cloud_index ! Initialize output arrays nclds(:,:) = 0 icld(:,:,:) = 0 cldtop(:,:,:) = 0 cldbas(:,:,:) = 0 ! Initialize work arrays ctop_work(:,:,:) = 0 cbas_work(:,:,:) = 0 ic_work(:,:,:) = mod(icld_k(:,:,:),100) ! locate distinct cloud layers, tops and bases do k=kcbtop,lk ! first locate tops and increment cloud layer counter, nclds ! check for no clouds above any type of cloud layer where ( ic_work(:,:,k).ge.1 .and. ic_work(:,:,k-1).eq.0) nclds(:,:) = nclds(:,:) + 1 ctop_work(:,:,k) = k endwhere ! check for stratiform cloud layer above convective cloud layer ! anvil and super anvil (types 6,7) are stratiform type clouds where ( (ic_work(:,:,k).ge.4 .and. ic_work(:,:,k).le.5) .and. & (ic_work(:,:,k-1).ge.1 .and. ic_work(:,:,k-1).le.3) .or. & (ic_work(:,:,k-1).ge.6) ) nclds(:,:) = nclds(:,:) + 1 ctop_work(:,:,k) = k endwhere ! check for convective cloud layer above stratiform cloud layer where ( (ic_work(:,:,k-1).ge.4 .and. ic_work(:,:,k-1).le.5) .and. & (ic_work(:,:,k).ge.1 .and. ic_work(:,:,k).le.3) .or. & (ic_work(:,:,k).ge.6) ) nclds(:,:) = nclds(:,:) + 1 ctop_work(:,:,k) = k endwhere if (k.lt.lk) then ! locate bases ! check for no clouds above any type of cloud layer where ( ic_work(:,:,k).ge.1 .and. ic_work(:,:,k+1).eq.0) cbas_work(:,:,k) = k endwhere ! check for stratiform cloud layer above convective cloud layer ! anvil and super anvil (types 6,7) are stratiform type clouds where ( (ic_work(:,:,k+1).ge.4 .and. ic_work(:,:,k+1).le.5) .and. & (ic_work(:,:,k).ge.1 .and. ic_work(:,:,k).le.3) .or. & (ic_work(:,:,k).ge.6) ) cbas_work(:,:,k) = k endwhere ! check for convective cloud layer above stratiform cloud layer where ( (ic_work(:,:,k).ge.4 .and. ic_work(:,:,k).le.5) .and. & (ic_work(:,:,k+1).ge.1 .and. ic_work(:,:,k+1).le.3) .or. & (ic_work(:,:,k+1).ge.6) ) cbas_work(:,:,k) = k endwhere endif do j = 1,jdim do i = 1,idim if (k.eq.lk .and. ic_work(i,j,k).ge.1) then cbas_work(i,j,k) = lk endif end do end do end do !----------------------------------------------------------------------- ! compress output cloud tops, bases and cloud markers to cloud levels do j = 1,jdim do i = 1,idim if (nclds(i,j) .ge.1) then kpmax = nclds(i,j) kp = 1 do k = kcbtop,lk if (ctop_work(i,j,k).ne. 0) then cldtop(i,j,kp) = ctop_work(i,j,k) icld(i,j,kp) = icld_k(i,j,k) endif if (cbas_work(i,j,k).ne. 0) then cldbas(i,j,kp) = cbas_work(i,j,k) kp = kp + 1 endif if (kp .gt. kpmax) go to 999 end do 999 continue endif end do end do !----------------------------------------------------------------------- ! error check for too many cloud layers ! find maximum number of cloud layers maxcld = maxval(nclds(:,:)) if (maxcld .gt. kx/2) then print *,'pe, NCLDS =', mpp_pe(),nclds call error_mesg ('diag_cloud, layr_top_base', & 'NCLDS too large', FATAL) endif !----------------------------------------------------------------------- end subroutine LAYR_TOP_BASE !####################################################################### subroutine CLDAMT_MN (nclds,cldtop,cldbas,phalf,camt,cldamt) !------------------------------------------------------------------- ! This subroutine computes mass-weighted vertical mean cloud amount ! for each distinct cloud layer !------------------------------------------------------------------- !=================================================================== ! Arguments (intent in) ! Namelist variables used in this routine (defined at top of module) ! NOFOG - logical switch for not allowing rhum clouds (or fog) ! to occur beneath a certain level (low_lev_cloud_index) -> ! nofog = true ! to allow clouds at the lowest model level -> nofog = false ! LOW_LEV_CLOUD_INDEX - level below which no clouds are allowed to occur ! when nofog = true (model level index) !------------------------------------------------------------------- integer, intent(out), dimension(:,:,:) :: cldtop,cldbas integer, intent(out), dimension(:,:) :: nclds real, intent(in), dimension (:,:,:) :: phalf ! NCLDS number of (random overlapping) clouds in column and also ! the current # for clouds to be operating on ! (dimensioned IDIM x JDIM ) ! CLDTOP index of cloud tops (at cloud levels) ! (dimensioned IDIM x JDIM x kx) ! CLDBAS index of cloud bottoms (at cloud levels) ! (dimensioned IDIM x JDIM x kx) ! PHALF Pressure at half model levels ! (dimensioned IDIM x JDIM x kx+1) !=================================================================== ! Arguments (intent inout) real, intent(inout), dimension (:,:,:) :: camt ! CAMT tentative merged cloud amounts ! (dimensioned IDIM x JDIM x kx) !======================================================================= real, intent(out), dimension(:,:,:) :: cldamt ! CLDAMT cloud amount (fraction) (at cloud levels) ! (dimensioned IDIM x JDIM x kx) !======================================================================= ! (Intent local) real, dimension(SIZE(camt,1),SIZE(camt,2),SIZE(camt,3)) :: delp,c_work,weight ! DELP pressure thickness of model layers ! (dimensioned IDIM x JDIM x kx) ! C_WORK cloud amount work array (at cloud levels) ! (dimensioned IDIM x JDIM x kx) ! WEIGHT vertical sum of press thickness over model levels for each ! cloud level (dimensioned IDIM x JDIM x kx) ! loop limits and indices integer idim,jdim,kx, kpr, lk, i, j, k integer maxcld !======================================================================= !======================================================================= idim = SIZE(camt,1) jdim = SIZE(camt,2) kx = SIZE(camt,3) ! define lower limitfor stratiform clouds if (nofog) then lk = low_lev_cloud_index else lk = kx endif ! find maximum number of cloud layers maxcld = maxval(nclds(:,:)) !----------------------------------------------------------------------- ! initialize output and work cloud amount arrays cldamt(:,:,:) = 0.0 c_work(:,:,:) = 0.0 weight(:,:,:) = 0.0 !----------------------------------------------------------------------- ! calculate pressure thickness of model layers, for vertcal averaging do k=1,kx delp (:,:,k) = phalf(:,:,k+1)-phalf(:,:,k) end do !----------------------------------------------------------------------- ! calculate mass weighted vertical mean cloud amount do kpr=1,maxcld do j=1,jdim do i=1,idim if (cldtop(i,j,kpr) .lt. cldbas(i,j,kpr)) then do k=cldtop(i,j,kpr),cldbas(i,j,kpr) c_work(i,j,kpr) = c_work(i,j,kpr) + & delp(i,j,k)*camt(i,j,k) weight(i,j,kpr) = weight(i,j,kpr) + delp(i,j,k) end do cldamt(i,j,kpr) = c_work(i,j,kpr)/weight(i,j,kpr) ! reset cloud amounts at model sigma levels to mass weighted for the layer do k=cldtop(i,j,kpr),cldbas(i,j,kpr) camt(i,j,k) = cldamt(i,j,kpr) end do else if ( (cldtop(i,j,kpr) .eq. cldbas(i,j,kpr)) .and. & (cldtop(i,j,kpr).gt.0 .and. cldtop(i,j,kpr).le.kx)) then k = cldtop(i,j,kpr) cldamt(i,j,kpr) = camt(i,j,k) else if ( (cldtop(i,j,kpr) .gt. cldbas(i,j,kpr)) .or. & (cldtop(i,j,kpr).lt.0) .or. (cldbas(i,j,kpr).lt.0)) then print *,'pe, i,j,kpr, cldtop,cldbas =', mpp_pe(),i,j,kpr, & cldtop(i,j,kpr),cldbas(i,j,kpr),cldtop(i,j,kpr),cldbas(i,j,kpr) call error_mesg ('diag_cloud, cldamt_mn', & 'invalid cldtop and/or cldbas', FATAL) endif end do end do end do end subroutine CLDAMT_MN !############################################################################# subroutine ANVIL_CIRRUS (lhighb,lmidb,nclds,cldtop,cldbas,cnvcntq,lgscldelq, & pfull,phalf,icld) !------------------------------------------------------------------------ ! This subroutine calculates anvil and super anvil cirrus cloud amounts !------------------------------------------------------------------------ !=================================================================== ! Arguments (intent in) !------------------------------------------------------------------- ! Namelist variables used in this routine (defined at top of module) ! HIGH_LEV_CLOUD_INDEX - level above which no clouds are allowed to form ! (model level index) !------------------------------------------------------------------- !------------------------------------------------------------------- ! parameters used in anvil_cirrus real, parameter :: cnvdelqmn = 0.0,lgscldelqmn = 1.0e-9 real, parameter :: lcldbas_prc = .55, lcldtop_prc = .225 ! CNVDELQMN - min threshold value of accumulated count of mixing ratio ! change due to convection above which convective clouds are ! allowed ! ( It is set = 0.0 as a parameter, because it not anticipated ! that it would be changed.) ! LGSCLDELQMN - min threshold value of averaged mixing ratio rate of ! change due to large scale condensation ! LCLDBAS_PRC - Threshold sigma value for base of an "anvil cirrus ! indicator" layer ! LCLDBAS_PRC - Threshold sigma value for top of an "anvil cirrus ! indicator" layer !------------------------------------------------------------------- integer, intent(out), dimension(:,:,:) :: cldtop,cldbas integer, intent(out), dimension(:,:) :: nclds,lhighb,lmidb real, intent(in), dimension (:,:,:) :: lgscldelq,cnvcntq,pfull,phalf !------------------------------------------------------------------- ! LHIGHB vertical level index lower limit for high cloud bases ! a function of lon and latitude (dimensioned IDIMxJDIM) ! LMIDB vertical level index lower limit for mid cloud bases ! a function of lon and latitude (dimensioned IDIMxJDIM) ! NCLDS number of (random overlapping) clouds in column and also ! the current # for clouds to be operating on ! (dimensioned IDIM x JDIM ) ! CLDTOP index of cloud tops (at cloud levels) ! (dimensioned IDIM x JDIM x kx) ! CLDBAS index of cloud bottoms (at cloud levels) ! (dimensioned IDIM x JDIM x kx) ! LGSCLDELQ Averaged rate of change in mix ratio due to lg scale precip ! at full model levels ! (dimensioned IDIM x JDIM x kx) ! CNVCNTQ Accumulated count of change in mix ratio due to conv precip ! at full model levels ! (dimensioned IDIM x JDIM x kx) ! PFULL Pressure at full model levels ! (dimensioned IDIM x JDIM x kx) ! PHALF Pressure at half model levels ! (dimensioned IDIM x JDIM x kx+1) !=================================================================== ! Arguments (intent inout) integer, intent(out), dimension (:,:,:) :: icld ! ICLD marker array of cloud types/heights (at cloud levels) ! (dimensioned IDIM x JDIM x kx) !======================================================================= ! (Intent local) integer, dimension (size(pfull,1),size(pfull,2),size(pfull,3)) :: ic_work !------------------------------------------------------------------- ! IC_WORK cloud marker work array ! (dimensioned IDIM x JDIM x kx) ! LCLDT Top of anvil cirrus ! (model level index) ! LCLDB Base of anvil cirrus ! (model level index) ! LPRCTOP Base of region in search for anvil cirrus markers ! (model level index) ! LPRCBAS Base of region in search for anvil cirrus markers ! (model level index) ! SIG_TOP_BASE Sigma value of top/base of region in search for ! anvil cirrus markers ! real sig_top_base integer lcldt, lcldb, lprcbas, lprctop integer i, j, idim, jdim, kx !----------------------------------------------------------------------- ! type loop index variable integer k,kl !=================================================================== idim = SIZE(pfull,1) jdim = SIZE(pfull,2) kx = size(pfull,3) ! Initialize cloud type work array ic_work (:,:,:) = icld(:,:,:) !----------------------------------------------------------------------- ! <><><><><><><><> search for anvil cirrus <><><><><><><><> !----------------------------------------------------------------------- do j=1,jdim do i=1,idim do k=1,kx sig_top_base = pfull(i,j,k)/phalf(i,j,kx+1) if ( sig_top_base .gt. lcldtop_prc) then lprctop = k-1 go to 100 endif end do 100 continue do k=1,kx sig_top_base = pfull(i,j,k)/phalf(i,j,kx+1) if ( sig_top_base .gt. lcldbas_prc) then lprcbas = k-1 go to 150 endif end do 150 continue do k=1,nclds(i,j) lcldt = cldtop(i,j,k) lcldb = cldbas(i,j,k) ! Must have a high cloud for any anvil cirrus if (icld(i,j,k) .ne. 101) go to 400 ! First check whether super-anvil cirrus criteria is satisfied ! This criteria differs from v197, because calculation of convective precip for ! high cloud region only is unreliable. Instead it is required that there has ! been a change in mixing ratio due to convection at at least one time step during ! the time interval over which the clouds are being calculated, at ALL levels ! between the base of the anvil cirrus indicator region (lprcbas- clculated from preset parameter) ! and the higher of either the actual top of the high cloud or the previoulsy computed top of the ! Anvil Cirrus cloud indicator region. if (lcldt .lt. lprctop) then lprctop = lcldt endif do kl=lprcbas,lprctop,-1 if (cnvcntq(i,j,kl) .lt. 1.0) then go to 200 endif end do ic_work(i,j,k) = 107 ! found super-anvil cirrus, no need to search for anvil cirrus go to 400 200 continue ! search for anvil cirrus using criteria that a change in mixing ratio due to ! either convection or large scale condensation at any level within an actual high ! cloud takes places during the time interval over which the clouds are being ! calculated. do kl=lcldb,lcldt,-1 if (cnvcntq(i,j,kl) .gt. 0.0) go to 300 if (kl.le.lhighb(i,j) .and. lgscldelq(i,j,kl).gt.lgscldelqmn) & go to 300 end do go to 400 300 continue ic_work(i,j,k) = 106 end do 400 continue ! endif end do end do !----------------------------------------------------------------------- ! <><><><><><><><> update cloud marker array <><><><><><><><> !----------------------------------------------------------------------- icld(:,:,:) = ic_work(:,:,:) !----------------------------------------------------------------------- end subroutine ANVIL_CIRRUS !####################################################################### subroutine CLD_LAYR_MN_TEMP_DELP & (nclds,cldtop,cldbas,temp,phalf,tempcld,delp_true) !------------------------------------------------------------------- ! This subroutine computes the mass-weighted mean cloud temperature ! and true cloud pressure thickness of distinct cloud layers. !------------------------------------------------------------------- !=================================================================== ! Arguments (intent in) ! Namelist variables used in this routine (defined at top of module) ! HIGH_LEV_CLOUD_INDEX - level above which no clouds are allowed to form ! (model level index) ! NOFOG - logical switch for not allowing rhum clouds (or fog) ! to occur beneath a certain level (low_lev_cloud_index) -> ! nofog = true ! to allow clouds at the lowest model level -> nofog = false ! LOW_LEV_CLOUD_INDEX - level below which no clouds are allowed to occur ! when nofog = true (model level index) !------------------------------------------------------------------- ! Arguments (intent in) integer, intent(in), dimension(:,:) :: nclds integer, intent(in), dimension(:,:,:) :: cldtop,cldbas real, intent(in), dimension(:,:,:) :: temp,phalf ! NCLDS number of (random overlapping) clouds in column and also ! the current # for clouds to be operating on ! (dimensioned IDIM x JDIM ) ! CLDTOP index of cloud tops (at cloud levels) ! (dimensioned IDIM x JDIM x kx) ! CLDBAS index of cloud bottoms (at cloud levels) ! (dimensioned IDIM x JDIM x kx) ! TEMP Temperature (Deg K) at full model levels ! (dimensioned IDIM x JDIM x kx) ! PHALF Pressure at half model levels ! (dimensioned IDIM x JDIM x kx+1) !=================================================================== ! Arguments (intent out) real, intent(out), dimension(:,:,:) :: tempcld,delp_true ! TEMPCLD cloud layer mean temperature (degrees Kelvin) ! (at cloud levels) ! DELP_TRUE true cloud pressure thickness of distinct cloud layers ! (at cloud levels) !======================================================================= ! (Intent local) integer, dimension (size(cldtop,1),size(cldtop,2)) :: cldt,cldb ! CLDT cloud top index work array ! (dimensioned IDIM x JDIM x kx) ! CLDB cloud base index work array ! (dimensioned IDIM x JDIM x kx) integer i,j,k,kk,kx,idim,jdim, lk, kcbtop, maxcld !======================================================================= !======================================================================= idim = SIZE(temp,1) jdim = SIZE(temp,2) kx = SIZE(temp,3) ! define lower limit of cloud bottom if (nofog) then lk = low_lev_cloud_index else lk = kx endif ! define upper limit of cloud tops kcbtop = high_lev_cloud_index ! Initialize output arrays tempcld(:,:,:) = 0 delp_true(:,:,:) = 0 ! find maximum number of clouds maxcld = maxval(nclds(:,:)) if (maxcld .ge. 1) then !----------------------------------------------------------------------- ! <><><><><><><><> calc pressure thickness <><><><><><><><> !----------------------------------------------------------------------- do k = 1,maxcld ! Initialize internal arrays cldt = 0 cldb = 0 where (nclds(:,:) .ge. k) cldt(:,:) = cldtop(:,:,k) cldb(:,:) = cldbas(:,:,k) end where do j=1,jdim do i=1,idim if (k .le. nclds(i,j)) then delp_true(i,j,k) = phalf(i,j,cldb(i,j)+1)-phalf(i,j,cldt(i,j)) endif end do end do end do !----------------------------------------------------------------------- ! <><><><><><><><> calc mass-weghted mean cloud temp <><><><><><><><> !----------------------------------------------------------------------- do j=1,jdim do i=1,idim do kk= 1,nclds(i,j) do k=cldtop(i,j,kk),cldbas(i,j,kk) tempcld(i,j,kk) = tempcld(i,j,kk) + & (phalf(i,j,k+1)-phalf(i,j,k))*temp(i,j,k)/delp_true(i,j,kk) end do end do end do end do endif !----------------------------------------------------------------------- end subroutine CLD_LAYR_MN_TEMP_DELP !####################################################################### SUBROUTINE DIAG_CLOUD_INIT( ix,iy,kx, ierr ) !======================================================================= ! ***** INITIALIZE Predicted Cloud Scheme !======================================================================= !--------------------------------------------------------------------- ! Arguments (Intent in) ! parmameter mxband = max number of radiative bands to be considered for some ! cloud properties (defined at top of module) !--------------------------------------------------------------------- integer, intent(in) :: ix, iy, kx ! INPUT ! ------ ! IX, IY, KX Dimensions for global storage arrays (2- horiz, vert) !--------------------------------------------------------------------- ! Arguments (Intent out) !--------------------------------------------------------------------- integer, intent(out) :: ierr ! OUTPUT ! ------ ! IERR Error flag !--------------------------------------------------------------------- ! (Intent local) !--------------------------------------------------------------------- integer unit, io, ierrnml, logunit integer id_restart character(len=32) :: fname !===================================================================== !--------------------------------------------------------------------- ! --- Read namelist !--------------------------------------------------------------------- if( FILE_EXIST( 'input.nml' ) ) then ! ------------------------------------- unit = open_namelist_file ('input.nml') ierrnml = 1 do while( ierrnml .ne. 0 ) READ ( unit, nml = diag_cloud_nml, iostat = io, end = 10 ) ierrnml = check_nml_error(io,'diag_cloud _nml') end do 10 call close_file (unit) ! ------------------------------------- end if !--------------------------------------------------------------------- ! --- Output namelist !--------------------------------------------------------------------- logunit = stdlog() if ( mpp_pe() == mpp_root_pe() ) then call write_version_number(version, tagname) write (logunit, nml=diag_cloud_nml) endif do_netcdf_restart = .true. call get_restart_io_mode(do_netcdf_restart) !--------------------------------------------------------------------- ! --- Allocate storage for global cloud quantities !--------------------------------------------------------------------- allocate( temp_sum(ix,iy,kx),qmix_sum(ix,iy,kx),rhum_sum(ix,iy,kx) ) allocate( qmix_sum2(ix,iy) ) allocate( omega_sum(ix,iy,kx),lgscldelq_sum(ix,iy,kx),cnvcntq_sum(ix,iy,kx) ) allocate( convprc_sum(ix,iy),nsum(ix,iy), nsum2(ix,iy) ) ! need to set up to account for first radiation step without having ! diag cloud info available (radiation called before diag_cloud, and ! diag_cloud being initiated (cold-started) in this job). tot_pts = ix*iy !--------------------------------------------------------------------- !---------- initialize for cloud averaging ------------------------- !--------------------------------------------------------------------- fname = 'diag_cloud.res.nc' if(do_netcdf_restart) then id_restart = register_restart_field(Dia_restart, fname, 'nsum', nsum, no_domain=.true.) id_restart = register_restart_field(Dia_restart, fname, 'temp_sum', temp_sum, no_domain=.true.) id_restart = register_restart_field(Dia_restart, fname, 'qmix_sum', qmix_sum, no_domain=.true.) id_restart = register_restart_field(Dia_restart, fname, 'rhum_sum', rhum_sum, no_domain=.true.) id_restart = register_restart_field(Dia_restart, fname, 'omega_sum', omega_sum, no_domain=.true.) id_restart = register_restart_field(Dia_restart, fname, 'lgscldelq_sum', lgscldelq_sum, no_domain=.true.) id_restart = register_restart_field(Dia_restart, fname, 'cnvcntq_sum', cnvcntq_sum, no_domain=.true.) id_restart = register_restart_field(Dia_restart, fname, 'convprc_sum', convprc_sum, no_domain=.true.) endif if( FILE_EXIST( 'INPUT/diag_cloud.res.nc' ) ) then if(mpp_pe() == mpp_root_pe() ) call error_mesg ('diag_cloud_mod', & 'Reading netCDF formatted restart file: INPUT/diag_cloud.res.nc', NOTE) call restore_state(Dia_restart) nsum2 = nsum qmix_sum2(:,:) = qmix_sum(:,:,size(qmix_sum,3)) ierr = 0 num_pts = tot_pts else if( FILE_EXIST( 'INPUT/diag_cloud.res' ) ) then unit = open_restart_file ('INPUT/diag_cloud.res', action='read') call read_data (unit,nsum) nsum2 = nsum call read_data (unit,temp_sum) call read_data (unit,qmix_sum) qmix_sum2(:,:) = qmix_sum(:,:,size(qmix_sum,3)) call read_data (unit,rhum_sum) call read_data (unit,omega_sum) call read_data (unit,lgscldelq_sum) call read_data (unit,cnvcntq_sum) call read_data (unit,convprc_sum) ierr = 0 num_pts = tot_pts else ierr = 1 if (mpp_pe() == mpp_root_pe() ) write (logunit,12) 12 format ('*** WARNING *** No cloud_tg restart file found *** ' ) nsum = 0 nsum2 = 0 temp_sum = 0.0 qmix_sum = 0.0 qmix_sum2 = 0.0 rhum_sum = 0.0 omega_sum = 0.0 lgscldelq_sum = 0.0 cnvcntq_sum = 0.0 convprc_sum = 0.0 num_pts = 0 end if !------------------------------------------------------------------- ! initialize zonal cloud routine for climatological zonal mean cloud info ! Passing the number 5 as the argument to cloud_zonal_init initializes the ! clim zonal clouds allowing seasonal variation call cloud_zonal_init (5) ! initialize shallow convection for shallow convective clouds call shallow_conv_init( kx ) ! Initialize cloud optical and radiative properties scheme if (.not. do_crad_init) then call diag_cloud_rad_init (do_crad_init) endif do_cpred_init = .true. module_is_initialized = .true. !===================================================================== end SUBROUTINE DIAG_CLOUD_INIT !####################################################################### SUBROUTINE DIAG_CLOUD_END !======================================================================= ! local integer :: unit !======================================================================= if( do_netcdf_restart) then if (mpp_pe() == mpp_root_pe()) then call error_mesg ('diag_cloud_mod', 'Writing netCDF formatted restart file: RESTART/diag_cloud.res.nc', NOTE) endif call diag_cloud_restart else if (mpp_pe() == mpp_root_pe()) then call error_mesg ('diag_cloud_mod', 'Writing native formatted restart file.', NOTE) endif unit = open_restart_file ('RESTART/diag_cloud.res', action='write') call write_data (unit, nsum) call write_data (unit, temp_sum) call write_data (unit, qmix_sum) call write_data (unit, rhum_sum) call write_data (unit, omega_sum) call write_data (unit, lgscldelq_sum) call write_data (unit, cnvcntq_sum) call write_data (unit, convprc_sum) call close_file (unit) endif module_is_initialized = .false. !===================================================================== end SUBROUTINE DIAG_CLOUD_END !####################################################################### ! ! ! ! write out restart file. ! Arguments: ! timestamp (optional, intent(in)) : A character string that represents the model time, ! used for writing restart. timestamp will append to ! the any restart file name as a prefix. ! ! subroutine diag_cloud_restart(timestamp) character(len=*), intent(in), optional :: timestamp if( do_netcdf_restart) then call save_restart(Dia_restart, timestamp) else call error_mesg ('diag_cloud_mod', & 'Native intermediate restart files are not supported.', FATAL) endif end subroutine diag_cloud_restart ! NAME="diag_cloud_restart" !####################################################################### function do_diag_cloud ( ) result (answer) logical :: answer ! returns logical value for whether diag_cloud has been initialized ! presumably if initialized then diag_cloud will be used answer = do_cpred_init end function do_diag_cloud !####################################################################### SUBROUTINE DIAG_CLOUD_SUM (is,js, & temp,qmix,rhum,omega,lgscldelq,cnvcntq,convprc,kbot) !----------------------------------------------------------------------- integer, intent(in) :: is,js real, intent(in), dimension (:,:,:) :: temp,qmix,rhum,omega real, intent(in), dimension (:,:,:) :: lgscldelq,cnvcntq real, intent(in), dimension (:,:) :: convprc integer, intent(in), OPTIONAL, dimension(:,:) :: kbot ! INPUT ! ------ ! IS,JS starting i,j indices from the full horizontal grid ! TEMP Temperature (Deg K) at full model levels ! (dimensioned IDIM x JDIM x kx) ! QMIX Mixing Ratio at full model levels ! (dimensioned IDIM x JDIM x kx) ! RHUM Relative humidity fraction at full model levels ! (dimensioned IDIM x JDIM x kx) ! OMEGA Pressure vertical velocity at full model levels ! (dimensioned IDIM x JDIM x kx) ! LGSCLDELQ Averaged rate of change in mix ratio due to lg scale precip ! at full model levels ! (dimensioned IDIM x JDIM x kx) ! CNVCNTQ Accumulated count of change in mix ratio due to conv precip ! at full model levels ! (dimensioned IDIM x JDIM x kx) ! convprc Accumulated conv precip rate summed over all ! full model levels (mm/day ) ! (dimensioned IDIM x JDIM) ! KBOT OPTIONAL; lowest model level index array ! (dimensioned IDIM x JDIM) ! ******* kbot will be used to select only those qmix values that are really ! ******* needed (typically this will be the bottom level except for ! ******* step mountains !----------------------------------------------------------------------- integer :: ie, je ie = is + size(rhum,1) - 1 je = js + size(rhum,2) - 1 !--------- use time-averaged or instantaneous clouds ----------- if (do_average) then nsum(is:ie,js:je) = nsum(is:ie,js:je) + 1 nsum2(is:ie,js:je) = nsum2(is:ie,js:je) + 1 temp_sum(is:ie,js:je,:) = temp_sum(is:ie,js:je,:) + temp(:,:,:) qmix_sum(is:ie,js:je,:) = qmix_sum(is:ie,js:je,:) + qmix(:,:,:) qmix_sum2(is:ie,js:je) = qmix_sum2(is:ie,js:je) + qmix(:,:,size(qmix,3)) rhum_sum(is:ie,js:je,:) = rhum_sum(is:ie,js:je,:) + rhum(:,:,:) omega_sum(is:ie,js:je,:) = omega_sum(is:ie,js:je,:) + omega(:,:,:) lgscldelq_sum(is:ie,js:je,:) = lgscldelq_sum(is:ie,js:je,:) & + lgscldelq(:,:,:) cnvcntq_sum(is:ie,js:je,:) = cnvcntq_sum(is:ie,js:je,:) + cnvcntq(:,:,:) convprc_sum(is:ie,js:je) = convprc_sum(is:ie,js:je) + convprc(:,:) else nsum(is:ie,js:je) = 1 nsum2(is:ie,js:je) = 1 temp_sum(is:ie,js:je,:) = temp(:,:,:) qmix_sum(is:ie,js:je,:) = qmix(:,:,:) qmix_sum2(is:ie,js:je) = qmix(:,:,size(qmix,3)) rhum_sum(is:ie,js:je,:) = rhum(:,:,:) omega_sum(is:ie,js:je,:) = omega(:,:,:) lgscldelq_sum(is:ie,js:je,:) = lgscldelq(:,:,:) cnvcntq_sum(is:ie,js:je,:) = cnvcntq(:,:,:) convprc_sum(is:ie,js:je) = convprc(:,:) endif !----------------------------------------------------------------------- end SUBROUTINE DIAG_CLOUD_SUM !####################################################################### subroutine DIAG_CLOUD_AVG (is, js, temp,qmix,rhum,omega, & lgscldelq,cnvcntq,convprc, ierr) !----------------------------------------------------------------------- integer, intent(in) :: is, js real, intent(inout), dimension(:,:,:) :: temp,qmix,rhum,omega real, intent(inout), dimension(:,:,:) :: lgscldelq,cnvcntq real, intent(inout), dimension(:,:) :: convprc integer, intent(out) :: ierr !----------------------------------------------------------------------- integer ::ie, je, num, k !----------------------------------------------------------------------- if (size(rhum,3) .ne. size(rhum_sum,3)) call error_mesg ( & 'diag_cloud_avg in diag_cloud_mod', & 'input argument has the wrong size',2) ie = is + size(rhum,1) - 1 je = js + size(rhum,2) - 1 num = count(nsum(is:ie,js:je) == 0) if (num > 0) then ! ----- no average, return error flag ----- !!! call error_mesg ('diag_cloud_avg in diag_cloud_mod', & !!! 'dividing by a zero counter', 2) ierr = 1 else ! ----- compute average ----- do k = 1, size(rhum,3) temp(:,:,k) = temp_sum(is:ie,js:je,k) / float(nsum(is:ie,js:je)) qmix(:,:,k) = qmix_sum(is:ie,js:je,k) / float(nsum(is:ie,js:je)) rhum(:,:,k) = rhum_sum(is:ie,js:je,k) / float(nsum(is:ie,js:je)) omega(:,:,k) = omega_sum(is:ie,js:je,k) / float(nsum(is:ie,js:je)) lgscldelq(:,:,k) = lgscldelq_sum(is:ie,js:je,k) / & float(nsum(is:ie,js:je)) enddo convprc(:,:) = convprc_sum(is:ie,js:je) / & float(nsum(is:ie,js:je)) ! The convective delta qmix count should be a sum, so no average is taken cnvcntq(:,:,:) = cnvcntq_sum(is:ie,js:je,:) ierr = 0 endif nsum(is:ie,js:je) = 0 temp_sum(is:ie,js:je,:) = 0.0 qmix_sum(is:ie,js:je,:) = 0.0 rhum_sum(is:ie,js:je,:) = 0.0 omega_sum(is:ie,js:je,:) = 0.0 lgscldelq_sum(is:ie,js:je,:) = 0.0 cnvcntq_sum(is:ie,js:je,:) = 0.0 convprc_sum(is:ie,js:je) = 0.0 !----------------------------------------------------------------------- end SUBROUTINE DIAG_CLOUD_AVG !####################################################################### subroutine DIAG_CLOUD_AVG2 (is, js, qmix, ierr) !----------------------------------------------------------------------- integer, intent(in) :: is, js real, intent(inout), dimension(:,:) :: qmix integer, intent(out) :: ierr !----------------------------------------------------------------------- integer ::ie, je, num !----------------------------------------------------------------------- ! if (size(qmix,3) .ne. size(qmix_sum2,3)) call error_mesg ( & ! 'diag_cloud_avg in diag_cloud_mod', & ! 'input argument has the wrong size',2) ie = is + size(qmix,1) - 1 je = js + size(qmix,2) - 1 num = count(nsum2(is:ie,js:je) == 0) if (num > 0) then ! ----- no average, return error flag ----- !!! call error_mesg ('diag_cloud_avg in diag_cloud_mod', & !!! 'dividing by a zero counter', 2) ierr = 1 else ! ----- compute average ----- qmix(:,:) = qmix_sum2(is:ie,js:je) / float(nsum2(is:ie,js:je)) ierr = 0 endif nsum2(is:ie,js:je) = 0 qmix_sum2(is:ie,js:je) = 0.0 !----------------------------------------------------------------------- end SUBROUTINE DIAG_CLOUD_AVG2 !####################################################################### end MODULE DIAG_CLOUD_MOD