diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml
index 6378596797..0011c2e942 100644
--- a/src/core_atmosphere/Registry.xml
+++ b/src/core_atmosphere/Registry.xml
@@ -784,6 +784,8 @@
+
+
@@ -934,10 +936,12 @@
+
+
@@ -2179,6 +2183,16 @@
description="configuration for convection schemes"
possible_values="`suite',`cu_kain_fritsch',`cu_tiedtke',`cu_ntiedtke',`cu_grell_freitas',`off'"/>
+
+
+
+
+
+
+
+
+
+
+
+
-
-
-
-
+
+
+
+
@@ -3100,7 +3130,10 @@
+ description="grid-scale horizontal cloud fraction"/>
+
+
@@ -3322,6 +3355,10 @@
description="tendency of cloud ice mixing ratio due to cumulus convection"
packages="cu_grell_freitas_in;cu_kain_fritsch_in;cu_ntiedtke_in"/>
+
+
-
-
diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_cloudiness.F b/src/core_atmosphere/physics/mpas_atmphys_driver_cloudiness.F
index 93e833b328..a85aa07606 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_driver_cloudiness.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_driver_cloudiness.F
@@ -49,6 +49,8 @@ module mpas_atmphys_driver_cloudiness
! cloud formation, but changes to the cloud water and cloud ice mixing ratios only affect the long wave and
! short wave radiation codes.
! Laura D. Fowler (laura@ucar.edu) / 2016-07-05.
+! * implemented cloud fraction and condensates from Kain-Fritsch convection parameterization.
+! Jerold A. Herwehe (herwehe.jerry@epa.gov) / 2016-12-08.
! * since we removed the local variable radt_cld_scheme from mpas_atmphys_vars.F, now defines radt_cld_scheme
! as a pointer to config_radt_cld_scheme.
! Laura D. Fowler (laura@ucar.edu) / 2017-02-16.
@@ -67,6 +69,11 @@ subroutine allocate_cloudiness
if(.not.allocated(dx_p) ) allocate(dx_p(ims:ime,jms:jme) )
if(.not.allocated(xland_p) ) allocate(xland_p(ims:ime,jms:jme) )
if(.not.allocated(cldfrac_p)) allocate(cldfrac_p(ims:ime,kms:kme,jms:jme))
+ if(.not.allocated(cldfracrad_p))allocate(cldfracrad_p(ims:ime,kms:kme,jms:jme))
+ if(.not.allocated(cldfrac_dp_p))allocate(cldfrac_dp_p(ims:ime,kms:kme,jms:jme))
+ if(.not.allocated(cldfrac_sh_p))allocate(cldfrac_sh_p(ims:ime,kms:kme,jms:jme))
+ if(.not.allocated(qccu_p) ) allocate(qccu_p(ims:ime,kms:kme,jms:jme) )
+ if(.not.allocated(qicu_p) ) allocate(qicu_p(ims:ime,kms:kme,jms:jme) )
if(.not.allocated(qvrad_p) ) allocate(qvrad_p(ims:ime,kms:kme,jms:jme) )
if(.not.allocated(qcrad_p) ) allocate(qcrad_p(ims:ime,kms:kme,jms:jme) )
if(.not.allocated(qirad_p) ) allocate(qirad_p(ims:ime,kms:kme,jms:jme) )
@@ -81,6 +88,11 @@ subroutine deallocate_cloudiness
if(allocated(dx_p) ) deallocate(dx_p )
if(allocated(xland_p) ) deallocate(xland_p )
if(allocated(cldfrac_p)) deallocate(cldfrac_p)
+ if(allocated(cldfracrad_p))deallocate(cldfracrad_p)
+ if(allocated(cldfrac_dp_p))deallocate(cldfrac_dp_p )
+ if(allocated(cldfrac_sh_p))deallocate(cldfrac_sh_p )
+ if(allocated(qccu_p) ) deallocate(qccu_p )
+ if(allocated(qicu_p) ) deallocate(qicu_p )
if(allocated(qvrad_p) ) deallocate(qvrad_p )
if(allocated(qcrad_p) ) deallocate(qcrad_p )
if(allocated(qirad_p) ) deallocate(qirad_p )
@@ -105,14 +117,16 @@ subroutine cloudiness_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
!local variables and pointers:
integer:: i,j,k
-
+ real(kind=RKIND),dimension(:,:),pointer:: cldfrac_dp,cldfrac_sh,qc_cu,qi_cu
real(kind=RKIND),pointer:: len_disp
real(kind=RKIND),dimension(:),pointer:: areaCell,meshDensity
real(kind=RKIND),dimension(:),pointer:: xland
+ character(len=StrKIND),pointer:: convection_scheme
!-----------------------------------------------------------------------------------------------------------------
call mpas_pool_get_config(configs,'config_len_disp',len_disp)
+ call mpas_pool_get_config(configs,'config_convection_scheme',convection_scheme)
call mpas_pool_get_array(mesh,'areaCell' ,areaCell )
call mpas_pool_get_array(mesh,'meshDensity',meshDensity)
@@ -126,28 +140,61 @@ subroutine cloudiness_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
dx_p(i,j) = dx_p(i,j)*0.001
xland_p(i,j) = xland(i)
enddo
-
- !--- initialize the radiative cloud fraction and water vapor, cloud water, cloud ice,
- ! and snow mixing ratios:
- do k = kts,kte
- do i = its,ite
- qvrad_p(i,k,j) = qv_p(i,k,j)
- qcrad_p(i,k,j) = qc_p(i,k,j)
- qirad_p(i,k,j) = qi_p(i,k,j)
- qsrad_p(i,k,j) = qs_p(i,k,j)
- cldfrac_p(i,k,j) = 0._RKIND
- enddo
- enddo
enddo
+ !--- initialize the radiative cloud fraction and water vapor, cloud water, cloud ice,
+ ! and snow mixing ratios:
+
+ convection_select: select case(convection_scheme)
+
+ case("cu_kain_fritsch")
+ call mpas_pool_get_array(diag_physics,'cldfrac_dp',cldfrac_dp)
+ call mpas_pool_get_array(diag_physics,'cldfrac_sh',cldfrac_sh)
+ call mpas_pool_get_array(diag_physics,'qc_cu' ,qc_cu )
+ call mpas_pool_get_array(diag_physics,'qi_cu' ,qi_cu )
+
+ do j = jts,jte
+ do k = kts,kte
+ do i = its,ite
+ qvrad_p(i,k,j) = qv_p(i,k,j)
+ qcrad_p(i,k,j) = qc_p(i,k,j)
+ qirad_p(i,k,j) = qi_p(i,k,j)
+ qsrad_p(i,k,j) = qs_p(i,k,j)
+ cldfrac_p(i,k,j) = 0._RKIND
+ cldfracrad_p(i,k,j) = 0._RKIND
+ cldfrac_dp_p(i,k,j) = cldfrac_dp(k,i)
+ cldfrac_sh_p(i,k,j) = cldfrac_sh(k,i)
+ qccu_p(i,k,j) = qc_cu(k,i)
+ qicu_p(i,k,j) = qi_cu(k,i)
+ enddo
+ enddo
+ enddo
+
+ case default
+ do j = jts,jte
+ do k = kts,kte
+ do i = its,ite
+ qvrad_p(i,k,j) = qv_p(i,k,j)
+ qcrad_p(i,k,j) = qc_p(i,k,j)
+ qirad_p(i,k,j) = qi_p(i,k,j)
+ qsrad_p(i,k,j) = qs_p(i,k,j)
+ cldfrac_p(i,k,j) = 0._RKIND
+ cldfracrad_p(i,k,j) = 0._RKIND
+ enddo
+ enddo
+ enddo
+
+ end select convection_select
+
end subroutine cloudiness_from_MPAS
!=================================================================================================================
- subroutine cloudiness_to_MPAS(diag_physics,its,ite)
+ subroutine cloudiness_to_MPAS(configs,diag_physics,its,ite)
!=================================================================================================================
!input arguments:
integer,intent(in):: its,ite
+ type(mpas_pool_type),intent(in) :: configs
!inout arguments:
type(mpas_pool_type),intent(inout):: diag_physics
@@ -155,19 +202,69 @@ subroutine cloudiness_to_MPAS(diag_physics,its,ite)
!local variables and pointers:
integer:: i,j,k
- real(kind=RKIND),dimension(:,:),pointer:: cldfrac
+ real(kind=RKIND),dimension(:,:),pointer:: cldfrac,cldfracwcu
+ real(kind=RKIND),dimension(:),pointer:: cldfract,cldfracwcut
+
+ real(kind=RKIND):: cldji,cldwcuji
+ character(len=StrKIND),pointer:: convection_scheme
!-----------------------------------------------------------------------------------------------------------------
- call mpas_pool_get_array(diag_physics,'cldfrac',cldfrac)
+ call mpas_pool_get_config(configs,'config_convection_scheme',convection_scheme)
+
+ call mpas_pool_get_array(diag_physics,'cldfrac' ,cldfrac)
+ call mpas_pool_get_array(diag_physics,'cldfract',cldfract)
+
+ convection_select: select case(convection_scheme)
+
+ case("cu_kain_fritsch")
+ call mpas_pool_get_array(diag_physics,'cldfracwcu',cldfracwcu)
+ call mpas_pool_get_array(diag_physics,'cldfracwcut',cldfracwcut)
+
+ do j = jts,jte
+ do k = kts,kte
+ do i = its,ite
+ cldfrac(k,i) = cldfrac_p(i,k,j)
+ cldfracwcu(k,i) = cldfracrad_p(i,k,j)
+ enddo
+ enddo
+ enddo
+
+ !compute cloud diagnosis (random overlapping) by ZCX (from WRF)
+ do j = jts,jte
+ do i = its,ite
+ cldji = 1._RKIND
+ cldwcuji = 1._RKIND
+ do k = kte-1,kts,-1
+ cldji = cldji * (1.0 - cldfrac_p(i,k,j))
+ cldwcuji = cldwcuji * (1.0 - cldfracrad_p(i,k,j))
+ enddo
+ cldfract(i) = 1.0 - cldji
+ cldfracwcut(i) = 1.0 - cldwcuji
+ enddo
+ enddo
- do j = jts,jte
- do k = kts,kte
- do i = its,ite
- cldfrac(k,i) = cldfrac_p(i,k,j)
- enddo
- enddo
- enddo
+ case default
+ do j = jts,jte
+ do k = kts,kte
+ do i = its,ite
+ cldfrac(k,i) = cldfrac_p(i,k,j)
+ enddo
+ enddo
+ enddo
+
+ !compute cloud diagnosis (random overlapping) by ZCX (from WRF)
+ do j = jts,jte
+ do i = its,ite
+ cldji = 1._RKIND
+ do k = kte-1,kts,-1
+ cldji = cldji * (1.0 - cldfrac_p(i,k,j))
+ enddo
+ cldfract(i) = 1.0 - cldji
+ enddo
+ enddo
+
+ end select convection_select
end subroutine cloudiness_to_MPAS
@@ -187,14 +284,20 @@ subroutine driver_cloudiness(configs,mesh,diag_physics,sfc_input,its,ite)
!local variables and pointers:
character(len=StrKIND),pointer:: radt_cld_scheme
+ logical, pointer:: cu_radt_feedback
integer:: i,j,k
+ real(kind=RKIND):: cldfrac_cu
+ character(len=StrKIND),pointer:: convection_scheme
+
!-----------------------------------------------------------------------------------------------------------------
!call mpas_log_write('')
!call mpas_log_write('--- enter subroutine driver_cloudiness:')
- call mpas_pool_get_config(configs,'config_radt_cld_scheme',radt_cld_scheme)
+ call mpas_pool_get_config(configs,'config_radt_cld_scheme' ,radt_cld_scheme )
+ call mpas_pool_get_config(configs,'config_convection_scheme',convection_scheme)
+ call mpas_pool_get_config(configs,'config_cu_radt_feedback' ,cu_radt_feedback )
!copy MPAS arrays to local arrays:
call cloudiness_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
@@ -226,8 +329,41 @@ subroutine driver_cloudiness(configs,mesh,diag_physics,sfc_input,its,ite)
end select cld_fraction_select
+ convection_select: select case(convection_scheme)
+
+ case("cu_kain_fritsch")
+ !--- assign radiative cloud fraction:
+ do j = jts,jte
+ do k = kts,kte
+ do i = its,ite
+ cldfracrad_p(i,k,j) = cldfrac_p(i,k,j)
+ enddo
+ enddo
+ enddo
+
+ if(cu_radt_feedback) then
+ !adjust cloud fraction and condensates before radiation calls if Cu feedback is enabled by
+ !incorporating subgrid-scale convective clouds into resolved cloudiness (based on CESM):
+ do j = jts,jte
+ do k = kts,kte
+ do i = its,ite
+ cldfrac_cu=cldfrac_dp_p(i,k,j)+cldfrac_sh_p(i,k,j) !Cu cloud fraction
+ cldfracrad_p(i,k,j)=(1.-cldfrac_cu)*cldfracrad_p(i,k,j) !update resolved cloud fraction for Cu punch-through
+ cldfracrad_p(i,k,j)=cldfracrad_p(i,k,j)+cldfrac_cu !new total cloud fraction for radiation
+ cldfracrad_p(i,k,j)=amin1(1.0,cldfracrad_p(i,k,j))
+ qcrad_p(i,k,j) = qcrad_p(i,k,j)+(qccu_p(i,k,j)*cldfrac_cu) !new total cloud water for radiation
+ qirad_p(i,k,j) = qirad_p(i,k,j)+(qicu_p(i,k,j)*cldfrac_cu) !new total cloud ice for radiation
+ enddo
+ enddo
+ enddo
+ endif
+
+ case default
+
+ end select convection_select
+
!copy local arrays to MPAS grid:
- call cloudiness_to_MPAS(diag_physics,its,ite)
+ call cloudiness_to_MPAS(configs,diag_physics,its,ite)
!call mpas_log_write('--- end subroutine driver_cloudiness.')
diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_convection.F b/src/core_atmosphere/physics/mpas_atmphys_driver_convection.F
index 763776bf7f..89c405f643 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_driver_convection.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_driver_convection.F
@@ -88,6 +88,8 @@ module mpas_atmphys_driver_convection
! * for the kain_fritsch parameterization of convection, change the definition of dx_p to match that used in the
! Grell-Freitas and "new Tiedtke" parameterization.
! Laura D. Fowler (laura@ucar.edu) / 2016-10-18.
+! * modified to accommodate Kain-Fritsch trigger option, cloud fraction, and condensates.
+! Jerold A. Herwehe (herwehe.jerry@epa.gov) / 2016-12-06.
! * since we removed the local variable convection_scheme from mpas_atmphys_vars.F, now defines convection_scheme
! as a pointer to config_convection_scheme.
! Laura D. Fowler (laura@ucar.edu) / 2017-02-16.
@@ -176,14 +178,19 @@ subroutine allocate_convection(configs)
if(.not.allocated(qicu_p) ) allocate(qicu_p(ims:ime,kms:kme,jms:jme) )
case ("cu_kain_fritsch")
- if(.not.allocated(dx_p) ) allocate(dx_p(ims:ime,jms:jme) )
- if(.not.allocated(area_p) ) allocate(area_p(ims:ime,jms:jme) )
- if(.not.allocated(nca_p) ) allocate(nca_p(ims:ime,jms:jme) )
- if(.not.allocated(cubot_p) ) allocate(cubot_p(ims:ime,jms:jme) )
- if(.not.allocated(cutop_p) ) allocate(cutop_p(ims:ime,jms:jme) )
- if(.not.allocated(w0avg_p) ) allocate(w0avg_p(ims:ime,kms:kme,jms:jme) )
- if(.not.allocated(rqrcuten_p) ) allocate(rqrcuten_p(ims:ime,kms:kme,jms:jme))
- if(.not.allocated(rqscuten_p) ) allocate(rqscuten_p(ims:ime,kms:kme,jms:jme))
+ if(.not.allocated(dx_p) ) allocate(dx_p(ims:ime,jms:jme) )
+ if(.not.allocated(area_p) ) allocate(area_p(ims:ime,jms:jme) )
+ if(.not.allocated(nca_p) ) allocate(nca_p(ims:ime,jms:jme) )
+ if(.not.allocated(cubot_p) ) allocate(cubot_p(ims:ime,jms:jme) )
+ if(.not.allocated(cutop_p) ) allocate(cutop_p(ims:ime,jms:jme) )
+ if(.not.allocated(w0avg_p) ) allocate(w0avg_p(ims:ime,kms:kme,jms:jme) )
+ if(.not.allocated(cldfrac_dp_p)) allocate(cldfrac_dp_p(ims:ime,kms:kme,jms:jme))
+ if(.not.allocated(cldfrac_sh_p)) allocate(cldfrac_sh_p(ims:ime,kms:kme,jms:jme))
+ if(.not.allocated(qccu_p) ) allocate(qccu_p(ims:ime,kms:kme,jms:jme) )
+ if(.not.allocated(qicu_p) ) allocate(qicu_p(ims:ime,kms:kme,jms:jme) )
+ if(.not.allocated(rqrcuten_p) ) allocate(rqrcuten_p(ims:ime,kms:kme,jms:jme) )
+ if(.not.allocated(rqscuten_p) ) allocate(rqscuten_p(ims:ime,kms:kme,jms:jme) )
+ if(.not.allocated(rqvdynten_p) ) allocate(rqvdynten_p(ims:ime,kms:kme,jms:jme) )
do i = its,ite
do j = jts,jte
@@ -192,11 +199,23 @@ subroutine allocate_convection(configs)
enddo
enddo
+ do i = its,ite
+ do k = kts,kte
+ do j = jts,jte
+ cldfrac_dp_p(i,k,j) = 0._RKIND
+ cldfrac_sh_p(i,k,j) = 0._RKIND
+ qccu_p(i,k,j) = 0._RKIND
+ qicu_p(i,k,j) = 0._RKIND
+ enddo
+ enddo
+ enddo
+
do i = its,ite
do k = kts,kte
do j = jts,jte
rqrcuten_p(i,k,j) = 0._RKIND
rqscuten_p(i,k,j) = 0._RKIND
+ rqvdynten_p(i,k,j)= 0._RKIND
enddo
enddo
enddo
@@ -326,8 +345,13 @@ subroutine deallocate_convection(configs)
if(allocated(cubot_p) ) deallocate(cubot_p )
if(allocated(cutop_p) ) deallocate(cutop_p )
if(allocated(w0avg_p) ) deallocate(w0avg_p )
+ if(allocated(cldfrac_dp_p) ) deallocate(cldfrac_dp_p )
+ if(allocated(cldfrac_sh_p) ) deallocate(cldfrac_sh_p )
+ if(allocated(qccu_p) ) deallocate(qccu_p )
+ if(allocated(qicu_p) ) deallocate(qicu_p )
if(allocated(rqrcuten_p) ) deallocate(rqrcuten_p )
if(allocated(rqscuten_p) ) deallocate(rqscuten_p )
+ if(allocated(rqvdynten_p) ) deallocate(rqvdynten_p )
case ("cu_tiedtke","cu_ntiedtke")
if(allocated(hfx_p) ) deallocate(hfx_p )
@@ -434,9 +458,9 @@ subroutine driver_convection(itimestep,configs,mesh,sfc_input,diag_physics,tend_
!variables specific to Kain_Fritsch parameterization:
logical:: warm_rain,adapt_step_flag
integer:: ktau
- real(kind=RKIND):: curr_secs
real(kind=RKIND):: cudt
real(kind=RKIND):: cudtacttime
+ integer, pointer:: trigger_kf
!CCPP-compliant flags:
character(len=StrKIND):: errmsg
@@ -464,7 +488,6 @@ subroutine driver_convection(itimestep,configs,mesh,sfc_input,diag_physics,tend_
cudt = dt_cu/60.
!... call to convection schemes:
- curr_secs = -1
cudtacttime = -1
adapt_step_flag = .false.
do j = jts, jte
@@ -514,6 +537,7 @@ subroutine driver_convection(itimestep,configs,mesh,sfc_input,diag_physics,tend_
call mpas_timer_stop('Grell-Freitas')
case ("cu_kain_fritsch")
+ call mpas_pool_get_config(configs, 'config_kfeta_trigger', trigger_kf)
if(itimestep == 1) then
ktau = itimestep
else
@@ -525,7 +549,7 @@ subroutine driver_convection(itimestep,configs,mesh,sfc_input,diag_physics,tend_
pcps = pres_hyd_p , t = t_p , &
dt = dt_dyn , ktau = ktau , &
dxCell = dx_p , areaCell = area_p , &
- cudt = cudt , curr_secs = curr_secs , &
+ cudt = cudt , &
adapt_step_flag = adapt_step_flag , rho = rho_p , &
raincv = raincv_p , pratec = pratec_p , &
nca = nca_p , u = u_p , &
@@ -546,6 +570,9 @@ subroutine driver_convection(itimestep,configs,mesh,sfc_input,diag_physics,tend_
rthcuten = rthcuten_p , rqvcuten = rqvcuten_p , &
rqccuten = rqccuten_p , rqrcuten = rqrcuten_p , &
rqicuten = rqicuten_p , rqscuten = rqscuten_p , &
+ rqvften = rqvdynten_p , trigger = trigger_kf , &
+ cldfra_dp_kf = cldfrac_dp_p , cldfra_sh_kf=cldfrac_sh_p, &
+ qc_kf = qccu_p , qi_kf = qicu_p , &
ids = ids , ide = ide , jds = jds , jde = jde , kds = kds , kde = kde , &
ims = ims , ime = ime , jms = jms , jme = jme , kms = kds , kme = kme , &
its = its , ite = ite , jts = jts , jte = jte , kts = kts , kte = kte &
@@ -644,6 +671,7 @@ subroutine convection_from_MPAS(dt_dyn,configs,mesh,sfc_input,diag_physics,tend_
real(kind=RKIND),dimension(:,:),pointer:: zgrid
real(kind=RKIND),dimension(:,:),pointer:: qc_cu,qi_cu
real(kind=RKIND),dimension(:,:),pointer:: w0avg
+ real(kind=RKIND),dimension(:,:),pointer:: cldfrac_dp,cldfrac_sh
real(kind=RKIND),dimension(:,:),pointer:: rthcuten,rqvcuten,rqccuten,rqicuten,rqrcuten,rqscuten
real(kind=RKIND),dimension(:,:),pointer:: rthblten,rthdynten,rthratenlw,rthratensw
real(kind=RKIND),dimension(:,:),pointer:: rqvblten,rqvdynten,rucuten,rvcuten
@@ -755,13 +783,18 @@ subroutine convection_from_MPAS(dt_dyn,configs,mesh,sfc_input,diag_physics,tend_
call mpas_pool_get_array(mesh,'areaCell',areaCell)
call mpas_pool_get_array(mesh,'meshDensity',meshDensity)
- call mpas_pool_get_array(diag_physics,'cubot' ,cubot )
- call mpas_pool_get_array(diag_physics,'cutop' ,cutop )
- call mpas_pool_get_array(diag_physics,'nca' ,nca )
- call mpas_pool_get_array(diag_physics,'w0avg' ,w0avg )
+ call mpas_pool_get_array(diag_physics,'cubot' ,cubot )
+ call mpas_pool_get_array(diag_physics,'cutop' ,cutop )
+ call mpas_pool_get_array(diag_physics,'nca' ,nca )
+ call mpas_pool_get_array(diag_physics,'w0avg' ,w0avg )
+ call mpas_pool_get_array(diag_physics,'cldfrac_dp',cldfrac_dp)
+ call mpas_pool_get_array(diag_physics,'cldfrac_sh',cldfrac_sh)
+ call mpas_pool_get_array(diag_physics,'qc_cu' ,qc_cu )
+ call mpas_pool_get_array(diag_physics,'qi_cu' ,qi_cu )
- call mpas_pool_get_array(tend_physics,'rqrcuten',rqrcuten)
- call mpas_pool_get_array(tend_physics,'rqscuten',rqscuten)
+ call mpas_pool_get_array(tend_physics,'rqrcuten' ,rqrcuten )
+ call mpas_pool_get_array(tend_physics,'rqscuten' ,rqscuten )
+ call mpas_pool_get_array(tend_physics,'rqvdynten',rqvdynten )
do j = jts,jte
do i = its,ite
@@ -771,9 +804,17 @@ subroutine convection_from_MPAS(dt_dyn,configs,mesh,sfc_input,diag_physics,tend_
cubot_p(i,j) = cubot(i)
cutop_p(i,j) = cutop(i)
+ do k = kts,kte
+ cldfrac_dp_p(i,k,j) = cldfrac_dp(k,i)
+ cldfrac_sh_p(i,k,j) = cldfrac_sh(k,i)
+ qccu_p(i,k,j) = qc_cu(k,i)
+ qicu_p(i,k,j) = qi_cu(k,i)
+ enddo
+
do k = kts,kte
rqrcuten_p(i,k,j) = rqrcuten(k,i)
rqscuten_p(i,k,j) = rqscuten(k,i)
+ rqvdynten_p(i,k,j)= rqvdynten(k,i)
enddo
!decreases the characteristic time period that convection remains active. When nca_p
@@ -881,6 +922,7 @@ subroutine convection_to_MPAS(configs,diag_physics,tend_physics,its,ite)
real(kind=RKIND),dimension(:),pointer :: xmb_total,xmb_shallow
real(kind=RKIND),dimension(:,:),pointer:: qc_cu,qi_cu
real(kind=RKIND),dimension(:,:),pointer:: w0avg
+ real(kind=RKIND),dimension(:,:),pointer:: cldfrac_dp,cldfrac_sh
real(kind=RKIND),dimension(:,:),pointer:: rthcuten,rqvcuten,rqccuten,rqicuten,rqrcuten,rqscuten
real(kind=RKIND),dimension(:,:),pointer:: rucuten,rvcuten
@@ -947,13 +989,17 @@ subroutine convection_to_MPAS(configs,diag_physics,tend_physics,its,ite)
enddo
case ("cu_kain_fritsch")
- call mpas_pool_get_array(diag_physics,'cubot',cubot)
- call mpas_pool_get_array(diag_physics,'cutop',cutop)
- call mpas_pool_get_array(diag_physics,'nca' ,nca )
- call mpas_pool_get_array(diag_physics,'w0avg',w0avg)
-
- call mpas_pool_get_array(tend_physics,'rqrcuten',rqrcuten)
- call mpas_pool_get_array(tend_physics,'rqscuten',rqscuten)
+ call mpas_pool_get_array(diag_physics,'cubot' ,cubot )
+ call mpas_pool_get_array(diag_physics,'cutop' ,cutop )
+ call mpas_pool_get_array(diag_physics,'nca' ,nca )
+ call mpas_pool_get_array(diag_physics,'w0avg' ,w0avg )
+ call mpas_pool_get_array(diag_physics,'cldfrac_dp',cldfrac_dp)
+ call mpas_pool_get_array(diag_physics,'cldfrac_sh',cldfrac_sh)
+ call mpas_pool_get_array(diag_physics,'qc_cu' ,qc_cu )
+ call mpas_pool_get_array(diag_physics,'qi_cu' ,qi_cu )
+
+ call mpas_pool_get_array(tend_physics,'rqrcuten',rqrcuten )
+ call mpas_pool_get_array(tend_physics,'rqscuten',rqscuten )
do j = jts,jte
do i = its,ite
@@ -961,9 +1007,13 @@ subroutine convection_to_MPAS(configs,diag_physics,tend_physics,its,ite)
cutop(i) = cutop_p(i,j)
nca(i) = nca_p(i,j)
do k = kts, kte
- w0avg(k,i) = w0avg_p(i,k,j)
- rqrcuten(k,i) = rqrcuten_p(i,k,j)
- rqscuten(k,i) = rqscuten_p(i,k,j)
+ w0avg(k,i) = w0avg_p(i,k,j)
+ cldfrac_dp(k,i) = cldfrac_dp_p(i,k,j)
+ cldfrac_sh(k,i) = cldfrac_sh_p(i,k,j)
+ qc_cu(k,i) = qccu_p(i,k,j)
+ qi_cu(k,i) = qicu_p(i,k,j)
+ rqrcuten(k,i) = rqrcuten_p(i,k,j)
+ rqscuten(k,i) = rqscuten_p(i,k,j)
enddo
enddo
enddo
@@ -1002,6 +1052,9 @@ subroutine update_convection_step1(configs,diag_physics,tend_physics,its,ite)
!local pointers:
character(len=StrKIND),pointer:: convection_scheme
real(kind=RKIND),dimension(:),pointer :: nca,cubot,cutop,cuprec,raincv
+ real(kind=RKIND),dimension(:,:),pointer:: qc_cu,qi_cu
+ real(kind=RKIND),dimension(:,:),pointer:: w0avg
+ real(kind=RKIND),dimension(:,:),pointer:: cldfrac_dp,cldfrac_sh
real(kind=RKIND),dimension(:,:),pointer:: rthcuten,rqvcuten,rqccuten,rqicuten,rqrcuten,rqscuten
!local variables and arrays:
@@ -1019,6 +1072,11 @@ subroutine update_convection_step1(configs,diag_physics,tend_physics,its,ite)
call mpas_pool_get_array(diag_physics,'cutop' ,cutop )
call mpas_pool_get_array(diag_physics,'cuprec',cuprec)
call mpas_pool_get_array(diag_physics,'raincv',raincv)
+ call mpas_pool_get_array(diag_physics,'w0avg' ,w0avg )
+ call mpas_pool_get_array(diag_physics,'cldfrac_dp',cldfrac_dp)
+ call mpas_pool_get_array(diag_physics,'cldfrac_sh',cldfrac_sh)
+ call mpas_pool_get_array(diag_physics,'qc_cu' ,qc_cu )
+ call mpas_pool_get_array(diag_physics,'qi_cu' ,qi_cu )
call mpas_pool_get_array(tend_physics,'rthcuten',rthcuten)
call mpas_pool_get_array(tend_physics,'rqvcuten',rqvcuten)
@@ -1036,6 +1094,11 @@ subroutine update_convection_step1(configs,diag_physics,tend_physics,its,ite)
if(nca(i) .lt. 0.5*dt_dyn) then
do k = kts,kte
+ w0avg(k,i) = 0._RKIND
+ cldfrac_dp(k,i) = 0._RKIND
+ cldfrac_sh(k,i) = 0._RKIND
+ qc_cu(k,i) = 0._RKIND
+ qi_cu(k,i) = 0._RKIND
rthcuten(k,i) = 0._RKIND
rqvcuten(k,i) = 0._RKIND
rqccuten(k,i) = 0._RKIND
diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_lw.F b/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_lw.F
index d4d271e50d..8b7d30da0c 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_lw.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_lw.F
@@ -82,6 +82,8 @@ module mpas_atmphys_driver_radiation_lw
! * in the call to rrtmg_lwrad, substituted the variables qv_p, qc_p, qi_p, and qs_p with qvrad_p, qcrad_p,
! qirad_p, and qsrad_p initialized in subroutine cloudiness_from_MPAS.
! Laura D. Fowler (laura@ucar.edu) / 2016-07-09.
+! * implemented cloud fraction and condensates feedback from Kain-Fritsch convection parameterization.
+! Jerold A. Herwehe (herwehe.jerry@epa.gov) / 2016-12-08.
! * substituted "use mpas_atmphys_o3climatology" with "use module_ra_rrtmg_vinterp" since we moved subroutine
! vinterp_ozn to is own module in physics_wrf.
! laura D. Fowler (laura@ucar.edu) / 2017-01-27.
@@ -315,16 +317,17 @@ subroutine radiation_lw_from_MPAS(xtime_s,configs,mesh,state,time_lev,diag_physi
type(mpas_pool_type),intent(inout):: diag_physics
!local pointers:
- logical,pointer:: config_o3climatology
- logical,pointer:: config_microp_re
+ logical,pointer:: config_o3climatology,cu_radt_feedback
character(len=StrKIND),pointer:: radt_lw_scheme
character(len=StrKIND),pointer:: microp_scheme
+ character(len=StrKIND),pointer:: convection_scheme
+ logical,pointer:: config_microp_re
real(kind=RKIND),dimension(:),pointer :: latCell,lonCell
real(kind=RKIND),dimension(:),pointer :: skintemp,snow,xice,xland
real(kind=RKIND),dimension(:),pointer :: m_ps,pin
real(kind=RKIND),dimension(:),pointer :: sfc_albedo,sfc_emiss
- real(kind=RKIND),dimension(:,:),pointer :: cldfrac,m_hybi,o3clim,o3vmr
+ real(kind=RKIND),dimension(:,:),pointer :: cldfrac,cldfracwcu,m_hybi,o3clim,o3vmr
real(kind=RKIND),dimension(:,:),pointer :: re_cloud,re_ice,re_snow
real(kind=RKIND),dimension(:,:,:),pointer:: aerosols,ozmixm
@@ -335,10 +338,12 @@ subroutine radiation_lw_from_MPAS(xtime_s,configs,mesh,state,time_lev,diag_physi
!-----------------------------------------------------------------------------------------------------------------
- call mpas_pool_get_config(configs,'config_microp_re' ,config_microp_re )
- call mpas_pool_get_config(configs,'config_o3climatology' ,config_o3climatology)
- call mpas_pool_get_config(configs,'config_radt_lw_scheme',radt_lw_scheme )
- call mpas_pool_get_config(configs,'config_microp_scheme' ,microp_scheme )
+ call mpas_pool_get_config(configs,'config_o3climatology' ,config_o3climatology)
+ call mpas_pool_get_config(configs,'config_radt_lw_scheme' ,radt_lw_scheme )
+ call mpas_pool_get_config(configs,'config_microp_scheme' ,microp_scheme )
+ call mpas_pool_get_config(configs,'config_microp_re' ,config_microp_re )
+ call mpas_pool_get_config(configs,'config_convection_scheme',convection_scheme )
+ call mpas_pool_get_config(configs,'config_cu_radt_feedback' ,cu_radt_feedback )
call mpas_pool_get_array(mesh,'latCell',latCell)
call mpas_pool_get_array(mesh,'lonCell',lonCell)
@@ -355,7 +360,6 @@ subroutine radiation_lw_from_MPAS(xtime_s,configs,mesh,state,time_lev,diag_physi
call mpas_pool_get_array(diag_physics,'sfc_albedo',sfc_albedo)
call mpas_pool_get_array(diag_physics,'sfc_emiss' ,sfc_emiss )
- call mpas_pool_get_array(diag_physics,'cldfrac' ,cldfrac )
call mpas_pool_get_array(diag_physics,'o3clim' ,o3clim )
call mpas_pool_get_array(diag_physics,'o3vmr' ,o3vmr )
call mpas_pool_get_array(diag_physics,'m_hybi' ,m_hybi )
@@ -370,13 +374,42 @@ subroutine radiation_lw_from_MPAS(xtime_s,configs,mesh,state,time_lev,diag_physi
xland_p(i,j) = xland(i)
enddo
enddo
- do j = jts,jte
- do k = kts,kte
- do i = its,ite
- cldfrac_p(i,k,j) = cldfrac(k,i)
- enddo
- enddo
- enddo
+
+ convection_select: select case(convection_scheme)
+
+ case("cu_kain_fritsch")
+ !assign appropriate cloud fraction:
+ if(cu_radt_feedback) then
+ call mpas_pool_get_array(diag_physics,'cldfracwcu',cldfracwcu)
+ do j = jts,jte
+ do k = kts,kte
+ do i = its,ite
+ cldfracrad_p(i,k,j) = cldfracwcu(k,i)
+ enddo
+ enddo
+ enddo
+ else
+ call mpas_pool_get_array(diag_physics,'cldfrac',cldfrac)
+ do j = jts,jte
+ do k = kts,kte
+ do i = its,ite
+ cldfracrad_p(i,k,j) = cldfrac(k,i)
+ enddo
+ enddo
+ enddo
+ endif
+
+ case default
+ call mpas_pool_get_array(diag_physics,'cldfrac',cldfrac)
+ do j = jts,jte
+ do k = kts,kte
+ do i = its,ite
+ cldfracrad_p(i,k,j) = cldfrac(k,i)
+ enddo
+ enddo
+ enddo
+
+ end select convection_select
!initialization:
do j = jts,jte
@@ -850,7 +883,7 @@ subroutine driver_radiation_lw(xtime_s,configs,mesh,state,time_lev,diag_physics,
! qv3d = qv_p , qc3d = qc_p , qi3d = qi_p , &
! qs3d = qs_p , cldfra3d = cldfrac_p , tsk = tsk_p , &
qv3d = qvrad_p , qc3d = qcrad_p , qi3d = qirad_p , &
- qs3d = qsrad_p , cldfra3d = cldfrac_p , tsk = tsk_p , &
+ qs3d = qsrad_p , cldfra3d = cldfracrad_p, tsk = tsk_p , &
emiss = sfc_emiss_p , xland = xland_p , xice = xice_p , &
snow = snow_p , icloud = icloud , o3input = o3input , &
noznlevels = num_oznlevels , pin = pin_p , o3clim = o3clim_p , &
@@ -908,7 +941,7 @@ subroutine driver_radiation_lw(xtime_s,configs,mesh,state,time_lev,diag_physics,
qg3d = qg_p , f_qc = f_qc , &
f_qr = f_qr , f_qi = f_qi , &
f_qs = f_qs , f_ice_phy = f_ice , &
- f_rain_phy = f_rain , cldfra = cldfrac_p , &
+ f_rain_phy = f_rain , cldfra = cldfracrad_p , &
xland = xland_p , xice = xice_p , &
num_months = num_months , levsiz = num_oznlevels , &
pin0 = pin_p , ozmixm = ozmixm_p , &
diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_sw.F b/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_sw.F
index bc2be079ef..96d1099d13 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_sw.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_sw.F
@@ -79,6 +79,8 @@ module mpas_atmphys_driver_radiation_sw
! * in the call to rrtmg_swrad, substituted the variables qv_p, qc_p, qi_p, and qs_p with qvrad_p, qcrad_p,
! qirad_p, and qsrad_p initialized in subroutine cloudiness_from_MPAS.
! Laura D. Fowler (laura@ucar.edu) / 2016-07-09.
+! * implemented cloud fraction and condensates feedback from Kain-Fritsch convection parameterization.
+! Jerold A. Herwehe (herwehe.jerry@epa.gov) / 2016-12-08.
! * in subroutines radiation_sw_from_MPAS, revised the initialization of re_cloud, re_ice, re_snow, to
! handle the case when the cloud microphysics parameterization is turned off, i.e. config_microp_scheme='off'.
! Laura D. Fowler (laura@ucar.edu) / 2017-02-10.
@@ -373,10 +375,11 @@ subroutine radiation_sw_from_MPAS(configs,mesh,state,time_lev,diag_physics,atm_i
integer:: i,j,k,n
!local pointers:
- logical,pointer:: config_o3climatology
- logical,pointer:: config_microp_re
+ logical,pointer:: config_o3climatology,cu_radt_feedback
character(len=StrKIND),pointer:: radt_sw_scheme
character(len=StrKIND),pointer:: microp_scheme
+ character(len=StrKIND),pointer:: convection_scheme
+ logical,pointer:: config_microp_re
real(kind=RKIND),dimension(:),pointer :: latCell,lonCell
real(kind=RKIND),dimension(:),pointer :: skintemp,snow,xice,xland
@@ -384,17 +387,19 @@ subroutine radiation_sw_from_MPAS(configs,mesh,state,time_lev,diag_physics,atm_i
real(kind=RKIND),dimension(:),pointer :: sfc_albedo,sfc_emiss
real(kind=RKIND),dimension(:),pointer :: taod5502d
real(kind=RKIND),dimension(:,:),pointer :: zgrid
- real(kind=RKIND),dimension(:,:),pointer :: cldfrac,m_hybi,o3clim
+ real(kind=RKIND),dimension(:,:),pointer :: cldfrac,cldfracwcu,m_hybi,o3clim
real(kind=RKIND),dimension(:,:),pointer :: re_cloud,re_ice,re_snow
real(kind=RKIND),dimension(:,:),pointer :: taod5503d
real(kind=RKIND),dimension(:,:,:),pointer:: aerosols,ozmixm
!-----------------------------------------------------------------------------------------------------------------
- call mpas_pool_get_config(configs,'config_microp_re' ,config_microp_re )
- call mpas_pool_get_config(configs,'config_o3climatology' ,config_o3climatology)
- call mpas_pool_get_config(configs,'config_radt_sw_scheme',radt_sw_scheme )
- call mpas_pool_get_config(configs,'config_microp_scheme' ,microp_scheme )
+ call mpas_pool_get_config(configs,'config_o3climatology' ,config_o3climatology)
+ call mpas_pool_get_config(configs,'config_radt_sw_scheme' ,radt_sw_scheme )
+ call mpas_pool_get_config(configs,'config_microp_scheme' ,microp_scheme )
+ call mpas_pool_get_config(configs,'config_microp_re' ,config_microp_re )
+ call mpas_pool_get_config(configs,'config_convection_scheme',convection_scheme )
+ call mpas_pool_get_config(configs,'config_cu_radt_feedback' ,cu_radt_feedback )
call mpas_pool_get_array(mesh,'latCell',latCell)
call mpas_pool_get_array(mesh,'lonCell',lonCell)
@@ -411,7 +416,6 @@ subroutine radiation_sw_from_MPAS(configs,mesh,state,time_lev,diag_physics,atm_i
call mpas_pool_get_array(diag_physics,'sfc_albedo',sfc_albedo)
call mpas_pool_get_array(diag_physics,'sfc_emiss' ,sfc_emiss )
- call mpas_pool_get_array(diag_physics,'cldfrac' ,cldfrac )
call mpas_pool_get_array(diag_physics,'o3clim' ,o3clim )
call mpas_pool_get_array(diag_physics,'m_hybi' ,m_hybi )
call mpas_pool_get_array(diag_physics,'m_ps' ,m_ps )
@@ -428,13 +432,42 @@ subroutine radiation_sw_from_MPAS(configs,mesh,state,time_lev,diag_physics,atm_i
xland_p(i,j) = xland(i)
enddo
enddo
- do j = jts,jte
- do k = kts,kte
- do i = its,ite
- cldfrac_p(i,k,j) = cldfrac(k,i)
- enddo
- enddo
- enddo
+
+ convection_select: select case(convection_scheme)
+
+ case("cu_kain_fritsch")
+ !assign appropriate cloud fraction:
+ if(cu_radt_feedback) then
+ call mpas_pool_get_array(diag_physics,'cldfracwcu',cldfracwcu)
+ do j = jts,jte
+ do k = kts,kte
+ do i = its,ite
+ cldfracrad_p(i,k,j) = cldfracwcu(k,i)
+ enddo
+ enddo
+ enddo
+ else
+ call mpas_pool_get_array(diag_physics,'cldfrac',cldfrac)
+ do j = jts,jte
+ do k = kts,kte
+ do i = its,ite
+ cldfracrad_p(i,k,j) = cldfrac(k,i)
+ enddo
+ enddo
+ enddo
+ endif
+
+ case default
+ call mpas_pool_get_array(diag_physics,'cldfrac',cldfrac)
+ do j = jts,jte
+ do k = kts,kte
+ do i = its,ite
+ cldfracrad_p(i,k,j) = cldfrac(k,i)
+ enddo
+ enddo
+ enddo
+
+ end select convection_select
!initialization:
do j = jts,jte
@@ -886,7 +919,7 @@ subroutine driver_radiation_sw(itimestep,configs,mesh,state,time_lev,diag_physic
! qv3d = qv_p , qc3d = qc_p , qi3d = qi_p , &
! qs3d = qs_p , cldfra3d = cldfrac_p , tsk = tsk_p , &
qv3d = qvrad_p , qc3d = qcrad_p , qi3d = qirad_p , &
- qs3d = qsrad_p , cldfra3d = cldfrac_p , tsk = tsk_p , &
+ qs3d = qsrad_p , cldfra3d = cldfracrad_p , tsk = tsk_p , &
albedo = sfc_albedo_p , xland = xland_p , xice = xice_p , &
snow = snow_p , coszr = coszr_p , xtime = xtime_m , &
gmt = gmt , julday = julday , radt = radt , &
@@ -936,7 +969,7 @@ subroutine driver_radiation_sw(itimestep,configs,mesh,state,time_lev,diag_physic
qg3d = qg_p , f_qc = f_qc , &
f_qr = f_qr , f_qi = f_qi , &
f_qs = f_qs , f_ice_phy = f_ice , &
- f_rain_phy = f_rain , cldfra = cldfrac_p , &
+ f_rain_phy = f_rain , cldfra = cldfracrad_p , &
xland = xland_p , xice = xice_p , &
num_months = num_months , levsiz = num_oznlevels , &
pin0 = pin_p , ozmixm = ozmixm_p , &
diff --git a/src/core_atmosphere/physics/mpas_atmphys_vars.F b/src/core_atmosphere/physics/mpas_atmphys_vars.F
index 5485f8fef8..8dff364727 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_vars.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_vars.F
@@ -70,6 +70,8 @@ module mpas_atmphys_vars
! * moved the declarations of arrays delta_p,wstar_p,uoce_p,and voce_p since they are now used in both modules
! module_bl_ysu.F and module_bl_mynn.F.
! Laura D. Fowler (laura@ucar.edu) / 2016-10-27.
+! * added cloud fraction and condensate arrays from Kain-Fritsch convection parameterization.
+! Jerold A. Herwehe (herwehe.jerry@epa.gov) / 2016-12-05.
! * added the variable opt_thcnd (option to treat thermal conductivity in NoahLSM). added additional options and
! arrays to run the Noah LSM scheme from WRF version 3.9.0.
! Laura D. Fowler (laura@ucar.edu) / 2017-01-27.
@@ -228,7 +230,7 @@ module mpas_atmphys_vars
znu_hyd_p !(pres_hyd_p / P0) needed in the Tiedtke convection scheme [Pa]
!=================================================================================================================
-!... variables related to ozone climatlogy:
+!... variables related to ozone climatology:
!=================================================================================================================
real(kind=RKIND),dimension(:,:,:),allocatable:: &
@@ -315,15 +317,22 @@ module mpas_atmphys_vars
rthcuten_p, &!
rqvcuten_p, &!
rqccuten_p, &!
- rqicuten_p !
+ rqicuten_p, &!
+ rqvdynten_p !
-!... kain fritsch specific arrays:
+!... kain-fritsch specific arrays:
real(kind=RKIND),dimension(:,:),allocatable:: &
cubot_p, &!lowest convective level [-]
cutop_p, &!highest convective level [-]
nca_p !counter for cloud relaxation time [-]
real(kind=RKIND),dimension(:,:,:),allocatable:: &
w0avg_p !
+ real(kind=RKIND),dimension(:,:,:),allocatable:: &
+ cldfrac_dp_p, &!cloud fraction from deep convection [-]
+ cldfrac_sh_p, &!cloud fraction from shallow convection [-]
+ cldfracwcu_p !cloud fraction including Cu [-]
+ real(kind=RKIND),dimension(:,:),allocatable:: &
+ cldfracwcut_p !total cloud fraction including Cu [-]
real(kind=RKIND),dimension(:,:,:),allocatable:: &
rqrcuten_p, &!
@@ -352,16 +361,18 @@ module mpas_atmphys_vars
xmb_shallow_p !
real(kind=RKIND),dimension(:,:,:),allocatable:: &
- rthdynten_p, &!
- qccu_p, &!
- qicu_p !
+ rthdynten_p !
real(kind=RKIND),dimension(:,:,:),allocatable:: &
rthraten_p !
+!... grell and kain-fritsch specific arrays:
+ real(kind=RKIND),dimension(:,:,:),allocatable:: &
+ qccu_p, &!convective cloud water mixing ratio [kg/kg]
+ qicu_p !convective cloud ice mixing ratio [kg/kg]
+
!... grell and tiedkte specific arrays:
real(kind=RKIND),dimension(:,:,:),allocatable:: &
- rqvdynten_p, &!
rqvdynblten_p, &!
rthdynblten_p !
@@ -750,11 +761,15 @@ module mpas_atmphys_vars
real(kind=RKIND),dimension(:,:,:),allocatable:: &
cldfrac_p, &!cloud fraction [-]
+ cldfracrad_p, &!cloud fraction local to cloudiness and radiation [-]
qvrad_p, &!water vapor mixing ratio local to cloudiness and radiation [kg/kg]
qcrad_p, &!cloud liquid water mixing ratio local to cloudiness and radiation [kg/kg]
qirad_p, &!cloud ice mixing ratio local to cloudiness and radiation [kg/kg]
qsrad_p !snow mixing ratio local to cloudiness and radiation [kg/kg]
+ real(kind=RKIND),dimension(:,:),allocatable:: &
+ cldfract_p !total cloud fraction [-]
+
!=================================================================================================================
!.. variables and arrays related to land-surface parameterization:
!=================================================================================================================
diff --git a/src/core_atmosphere/physics/physics_wrf/module_cu_kfeta.F b/src/core_atmosphere/physics/physics_wrf/module_cu_kfeta.F
index ecf0d82adf..9f48c5e218 100644
--- a/src/core_atmosphere/physics/physics_wrf/module_cu_kfeta.F
+++ b/src/core_atmosphere/physics/physics_wrf/module_cu_kfeta.F
@@ -1,10 +1,47 @@
+!===============================================================================
+! This version of module_cu_kfeta.F comes from WRF v3.8.1 (also the same
+! module as in WRF v4.4) and was adapted for use in MPAS, including the KF
+! trigger option and subgrid-scale cloudiness and condensate feedbacks to the
+! radiation schemes.
+! Jerold A. Herwehe (herwehe.jerry@epa.gov) / 2016-12-05.
+!
+! This version of module_cu_kfeta.F also computes the convective time
+! scale (variable TIMEC) based on updraft dynamics, rather than a
+! horizontal grid scale advective time. Look-up table error checking
+! is also improved.
+! O. Russell Bullock Jr. (bullock.russell@epa.gov) / 2017-01-04.
+!
+! Additional modifications:
+! * Generalized module for use in both the MPAS and WRF models.
+! Jerold A. Herwehe (herwehe.jerry@epa.gov) / 2019-11-19.
+!===============================================================================
+!
MODULE module_cu_kfeta
-! USE module_wrf_error
+#if defined(mpas)
+use mpas_atmphys_utilities, only: physics_message,physics_error_fatal
+#define FATAL_ERROR(M) call physics_error_fatal( M )
+#define WRITE_MESSAGE(M) call physics_message( M )
+#else
+use module_wrf_error
+#define FATAL_ERROR(M) call wrf_error_fatal( M )
+#define WRITE_MESSAGE(M) call wrf_message( M )
+#endif
+!
+! V3.3: A new trigger function is added based Ma and Tan (2009):
+! Ma, L.-M. and Z.-M. Tan, 2009: Improving the behavior of
+! the cumulus parameterization for tropical cyclone prediction:
+! Convection trigger. Atmospheric Research, 92, 190 - 211.
+!
+!ckay
+! WRF v3.5 with diagnosed deep and shallow KF cloud fraction using
+! CAM3-CAM5 methodology, along with captured liquid and ice condensates.
+! JAH & KA (U.S. EPA) -- May 2013
+!
!--------------------------------------------------------------------
! Lookup table variables:
- INTEGER, PARAMETER :: KFNT=250,KFNP=220
+ INTEGER, PARAMETER, PRIVATE :: KFNT=300,KFNP=250
REAL, DIMENSION(KFNT,KFNP),PRIVATE, SAVE :: TTAB,QSTAB
REAL, DIMENSION(KFNP),PRIVATE, SAVE :: THE0K
REAL, DIMENSION(200),PRIVATE, SAVE :: ALU
@@ -16,38 +53,50 @@ MODULE module_cu_kfeta
CONTAINS
#if defined(mpas)
- SUBROUTINE KF_eta_CPS( &
- ids,ide, jds,jde, kds,kde &
- ,ims,ime, jms,jme, kms,kme &
- ,its,ite, jts,jte, kts,kte &
- ,DT,KTAU,dxCell,areaCell,CUDT,CURR_SECS,ADAPT_STEP_FLAG &
- ,rho,RAINCV,PRATEC,NCA &
- ,U,V,TH,T,W,dz8w,Pcps,pi &
- ,W0AVG,XLV0,XLV1,XLS0,XLS1,CP,R,G,EP1 &
- ,EP2,SVP1,SVP2,SVP3,SVPT0 &
- ,STEPCU,CU_ACT_FLAG,warm_rain,CUTOP,CUBOT &
- ,QV &
+ SUBROUTINE KF_eta_CPS( &
+ ids,ide, jds,jde, kds,kde &
+ ,ims,ime, jms,jme, kms,kme &
+ ,its,ite, jts,jte, kts,kte &
+ ,trigger &
+ ,DT,KTAU,dxCell,areaCell,CUDT,ADAPT_STEP_FLAG &
+ ,rho,RAINCV,PRATEC,NCA &
+ ,U,V,TH,T,W,dz8w,Pcps,pi &
+ ,W0AVG,XLV0,XLV1,XLS0,XLS1,CP,R,G,EP1 &
+ ,EP2,SVP1,SVP2,SVP3,SVPT0 &
+ ,STEPCU,CU_ACT_FLAG,warm_rain,CUTOP,CUBOT &
+ ,QV &
! optionals
- ,F_QR ,F_QI ,F_QS &
- ,RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN &
- ,RQICUTEN,RQSCUTEN &
+ ,F_QR ,F_QI ,F_QS &
+ ,RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN &
+ ,RQICUTEN,RQSCUTEN, RQVFTEN &
+ ,cldfra_dp_KF,cldfra_sh_KF &
+ ,qc_KF,qi_KF &
)
#else
SUBROUTINE KF_eta_CPS( &
ids,ide, jds,jde, kds,kde &
,ims,ime, jms,jme, kms,kme &
,its,ite, jts,jte, kts,kte &
- ,DT,KTAU,DX,CUDT,CURR_SECS,ADAPT_STEP_FLAG &
+ ,trigger &
+ ,DT,KTAU,DX,CUDT,ADAPT_STEP_FLAG &
,rho,RAINCV,PRATEC,NCA &
,U,V,TH,T,W,dz8w,Pcps,pi &
,W0AVG,XLV0,XLV1,XLS0,XLS1,CP,R,G,EP1 &
,EP2,SVP1,SVP2,SVP3,SVPT0 &
,STEPCU,CU_ACT_FLAG,warm_rain,CUTOP,CUBOT &
,QV &
+ ,shall & !BSINGH - For WRFCuP scheme added "shall"
! optionals
- ,F_QR ,F_QI ,F_QS &
+ ,F_QV ,F_QC ,F_QR ,F_QI ,F_QS &
,RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN &
- ,RQICUTEN,RQSCUTEN &
+ ,RQICUTEN,RQSCUTEN, RQVFTEN &
+!ckay
+ ,cldfra_dp_KF,cldfra_sh_KF &
+ ,qc_KF,qi_KF &
+!kf_edrates
+ ,UDR_KF,DDR_KF &
+ ,UER_KF,DER_KF &
+ ,TIMEC_KF,KF_EDRATES &
)
#endif
@@ -60,6 +109,7 @@ SUBROUTINE KF_eta_CPS( &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte
+ INTEGER, INTENT(IN ) :: trigger
INTEGER, INTENT(IN ) :: STEPCU
LOGICAL, INTENT(IN ) :: warm_rain
@@ -84,7 +134,7 @@ SUBROUTINE KF_eta_CPS( &
!
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
INTENT(INOUT) :: &
- W0AVG
+ W0AVG
!MPAS specific (Laura D. Fowler 2013-08-23): In MPAS, the mean distance between cell
!centers varies between cells, in particular, for variable-resolution meshes. Below,
@@ -97,8 +147,7 @@ SUBROUTINE KF_eta_CPS( &
#endif
REAL, INTENT(IN ) :: CUDT
- REAL, INTENT(IN ) :: CURR_SECS
- LOGICAL,INTENT(IN ) :: ADAPT_STEP_FLAG
+ LOGICAL,OPTIONAL,INTENT(IN ) :: ADAPT_STEP_FLAG
!
REAL, DIMENSION( ims:ime , jms:jme ), &
INTENT(INOUT) :: RAINCV
@@ -109,6 +158,11 @@ SUBROUTINE KF_eta_CPS( &
REAL, DIMENSION( ims:ime , jms:jme ), &
INTENT(INOUT) :: NCA
+#if !(defined(mpas))
+ REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
+ INTENT(INOUT) :: SHALL !BSINGH - For WRFCuP scheme
+#endif
+
REAL, DIMENSION( ims:ime , jms:jme ), &
INTENT(OUT) :: CUBOT, &
CUTOP
@@ -128,8 +182,8 @@ SUBROUTINE KF_eta_CPS( &
RQCCUTEN, &
RQRCUTEN, &
RQICUTEN, &
- RQSCUTEN
-
+ RQSCUTEN, &
+ RQVFTEN
!
! Flags relating to the optional tendency arrays declared above
! Models that carry the optional tendencies will provdide the
@@ -142,6 +196,29 @@ SUBROUTINE KF_eta_CPS( &
,F_QI &
,F_QS
+!ckay
+ REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
+ INTENT(INOUT) :: &
+ cldfra_dp_KF, &
+ cldfra_sh_KF, &
+ qc_KF, &
+ qi_KF
+
+#if !(defined(mpas))
+!kf_edrates
+ REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
+ INTENT(INOUT) :: &
+ UDR_KF, &
+ DDR_KF, &
+ UER_KF, &
+ DER_KF
+
+ REAL, DIMENSION( ims:ime , jms:jme ) , &
+ INTENT(INOUT) :: &
+ TIMEC_KF
+
+ INTEGER, INTENT(IN) :: KF_EDRATES
+#endif
! LOCAL VARS
@@ -155,6 +232,8 @@ SUBROUTINE KF_eta_CPS( &
QV1D, &
P1D, &
RHO1D, &
+ tpart_v1D, &
+ tpart_h1D, &
W0AVG1D
REAL, DIMENSION( kts:kte ):: &
@@ -165,23 +244,36 @@ SUBROUTINE KF_eta_CPS( &
DQSDT, &
DTDT
+ REAL, DIMENSION (its-1:ite+1,kts:kte,jts-1:jte+1) :: aveh_t, aveh_q
+ REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: aveh_qmax, aveh_qmin
+ REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: avev_t, avev_q
+ REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: avev_qmax, avev_qmin
+ REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: coef_v, coef_h, tpart_h, tpart_v
+ INTEGER :: ii,jj,kk
+
+ REAL :: ttop
+ REAL, DIMENSION (kts:kte) :: z0
+
#if defined(mpas)
real:: tst,tv,prs,rhoe,w0,scr1,dx,dxsq,tmp
#else
REAL :: TST,tv,PRS,RHOE,W0,SCR1,DXSQ,tmp
#endif
+ integer :: ibegh,iendh,jbegh,jendh
+ integer :: istart,iend,jstart,jend
INTEGER :: i,j,k,NTST
+#if !(defined(mpas))
+ INTEGER :: ishall !BSINGH - For WRFCuP Scheme
+#endif
REAL :: lastdt = -1.0
REAL :: W0AVGfctr, W0fctr, W0den
- LOGICAL :: run_param
!
#if !(defined(mpas))
DXSQ=DX*DX
#endif
-
!----------------------
NTST=STEPCU
TST=float(NTST*2)
@@ -227,33 +319,152 @@ SUBROUTINE KF_eta_CPS( &
ENDDO
lastdt = dt
-
-!
-!...CHECK FOR CONVECTIVE INITIATION EVERY 5 MINUTES (OR NTST/2)...
-!
-!----------------------
+! New trigger function
+ IF (trigger.eq.2) THEN
+!
+! calculate 9-point average of moisture advection and temperature using halo (Horizontal)
+!
+ aveh_t=-999 ! horizontal 9-point ave
+ aveh_q=-999
+ avev_t=0 ! vertical 3-level ave
+ avev_q=0
+ avev_qmax=0
+ avev_qmin=0
+ aveh_qmax=0
+ aveh_qmin=0
+ tpart_h=0
+ tpart_v=0
+ coef_h=0
+ coef_v=0
+ ibegh=max(its-1, ids+1) ! start from 2
+ jbegh=max(jts-1, jds+1)
+ iendh=min(ite+1, ide-2) ! end at ide-2
+ jendh=min(jte+1, jde-2)
+ DO J = jbegh,jendh
+ DO K = kts,kte
+ DO I = ibegh,iendh
+ aveh_t(i,k,j)=(T(i-1,k,j-1)+T(i-1,k,j) +T(i-1,k,j+1)+ &
+ T(i,k,j-1) +T(i,k,j) +T(i,k,j+1)+ &
+ T(i+1,k,j-1) +T(i+1,k,j) +T(i+1,k,j+1))/9.
+ aveh_q(i,k,j)=(rqvften(i-1,k,j-1)+rqvften(i-1,k,j) +rqvften(i-1,k,j+1)+ &
+ rqvften(i,k,j-1) +rqvften(i,k,j) +rqvften(i,k,j+1)+ &
+ rqvften(i+1,k,j-1) +rqvften(i+1,k,j) +rqvften(i+1,k,j+1))/9.
+ ENDDO
+ ENDDO
+ ENDDO
+! boundary value ( all processors will do the following? Or just those processsors handling sub-area including boundary)
+ DO K = kts,kte
+ DO J = jts-1,jte+1
+ DO I = its-1,ite+1
+
+ if(i.eq.ids) then
+ aveh_t(i,k,j)=aveh_t(i+1,k,j)
+ aveh_q(i,k,j)=aveh_q(i+1,k,j)
+ elseif(i.eq.ide-1) then
+ aveh_t(i,k,j)=aveh_t(i-1,k,j)
+ aveh_q(i,k,j)=aveh_q(i-1,k,j)
+ endif
-!
-! Modified for adaptive time step
-!
- if (ADAPT_STEP_FLAG) then
- if ( (KTAU .eq. 1) .or. (cudt .eq. 0) .or. &
- ( CURR_SECS + dt >= &
- ( int( CURR_SECS / ( cudt * 60 ) ) + 1 ) * cudt * 60 ) ) then
- run_param = .TRUE.
- else
- run_param = .FALSE.
- endif
+ if(j.eq.jds) then
+ aveh_t(i,k,j)=aveh_t(i,k,j+1)
+ aveh_q(i,k,j)=aveh_q(i,k,j+1)
+ elseif(j.eq.jde-1) then
+ aveh_t(i,k,j)=aveh_t(i,k,j-1)
+ aveh_q(i,k,j)=aveh_q(i,k,j-1)
+ endif
- else
- if (MOD(KTAU,NTST) .EQ. 0 .or. KTAU .eq. 1) then
- run_param = .TRUE.
- else
- run_param = .FALSE.
- endif
- endif
+ if(j.eq.jds.and.i.eq.ids) then
+ aveh_q(i,k,j)=aveh_q(i+1,k,j+1)
+ aveh_t(i,k,j)=aveh_t(i+1,k,j+1)
+ endif
- if (run_param) then
+ if(j.eq.jde-1.and.i.eq.ids) then
+ aveh_q(i,k,j)=aveh_q(i+1,k,j-1)
+ aveh_t(i,k,j)=aveh_t(i+1,k,j-1)
+ endif
+
+ if(j.eq.jde-1.and.i.eq.ide-1) then
+ aveh_q(i,k,j)=aveh_q(i-1,k,j-1)
+ aveh_t(i,k,j)=aveh_t(i-1,k,j-1)
+ endif
+
+ if(j.eq.jds.and.i.eq.ide-1) then
+ aveh_q(i,k,j)=aveh_q(i-1,k,j+1)
+ aveh_t(i,k,j)=aveh_t(i-1,k,j+1)
+ endif
+
+ ENDDO
+ ENDDO
+ ENDDO
+! search for max/min moisture advection in 9-point range, calculate horizontal T-perturbation (tpart_h)
+ istart=max(its, ids+1) ! start from 2
+ jstart=max(jts, jds+1)
+ iend=min(ite, ide-2) ! end at ide-2
+ jend=min(jte, jde-2)
+ DO K = kts,kte
+ DO J = jstart,jend
+ DO I = istart,iend
+ aveh_qmax(i,k,j)=aveh_q(i,k,j)
+ aveh_qmin(i,k,j)=aveh_q(i,k,j)
+ DO ii=-1, 1
+ DO jj=-1,1
+ if(aveh_q(i+II,k,j+JJ).gt.aveh_qmax(i,k,j)) aveh_qmax(i,k,j)=aveh_q(i+II,k,j+JJ)
+ if(aveh_q(i+II,k,j+JJ).lt.aveh_qmin(i,k,j)) aveh_qmin(i,k,j)=aveh_q(i+II,k,j+JJ)
+ ENDDO
+ ENDDO
+ if(aveh_qmax(i,k,j).gt.aveh_qmin(i,k,j))then
+ coef_h(i,k,j)=(aveh_q(i,k,j)-aveh_qmin(i,k,j))/(aveh_qmax(i,k,j)-aveh_qmin(i,k,j))
+ else
+ coef_h(i,k,j)=0.
+ endif
+ coef_h(i,k,j)=amin1(coef_h(i,k,j),1.0)
+ coef_h(i,k,j)=amax1(coef_h(i,k,j),0.0)
+ tpart_h(i,k,j)=coef_h(i,k,j)*(T(i,k,j)-aveh_t(i,k,j))
+ ENDDO
+ ENDDO
+ ENDDO
+ 89 continue
+! vertical 3-layer calculation
+ DO J = jts, jte
+ DO I = its, ite
+ z0(1) = 0.5 * dz8w(i,1,j)
+ DO K = 2, kte
+ Z0(K) = Z0(K-1) + .5 * (DZ8W(i,K,j) + DZ8W(i,K-1,j))
+ ENDDO
+ DO K = kts+1,kte-1
+ ttop = t(i,k,j) + ((t(i,k,j) - t(i,k+1,j)) / (z0(k) - z0(k+1))) * (z0(k)-z0(k-1))
+ avev_t(i,k,j)=(T(i,k-1,j) + T(i,k,j) + ttop)/3.
+! avev_t(i,k,j)=(T(i,k-1,j)+T(i,k,j) + T(i,k+1,j))/3.
+ avev_q(i,k,j)=(rqvften(i,k-1,j)+rqvften(i,k,j) + rqvften(i,k+1,j))/3.
+ ENDDO
+ avev_t(i,kts,j)=avev_t(i,kts+1,j) ! lowest level value, is it the same as avev_t(i,kds,j)=avev_t(i,kds+1,j)?
+ avev_q(i,kts,j)=avev_q(i,kts+1,j)
+ avev_t(i,kte,j)=avev_t(i,kte-1,j) ! highest level value
+ avev_q(i,kte,j)=avev_q(i,kte-1,j)
+ ENDDO
+ ENDDO
+! max /min value
+ DO J = jts, jte
+ DO I = its, ite
+ DO K = kts+1,kte-1
+ avev_qmax(i,k,j)=avev_q(i,k,j)
+ avev_qmin(i,k,j)=avev_q(i,k,j)
+ DO kk=-1,1
+ if(avev_q(i,k+kk,j).gt.avev_qmax(i,k,j)) avev_qmax(i,k,j)=avev_q(i,k+kk,j)
+ if(avev_q(i,k+kk,j).lt.avev_qmin(i,k,j)) avev_qmin(i,k,j)=avev_q(i,k+kk,j)
+ ENDDO
+ if(avev_qmax(i,k,j).gt.avev_qmin(i,k,j)) then
+ coef_v(i,k,j)=(avev_q(i,k,j)-avev_qmin(i,k,j))/(avev_qmax(i,k,j)-avev_qmin(i,k,j))
+ else
+ coef_v(i,k,j)=0
+ endif
+ tpart_v(i,k,j)=coef_v(i,k,j)*(T(i,k,j)-avev_t(i,k,j))
+ ENDDO
+ tpart_v(i,kts,j)= tpart_v(i,kts+1,j) ! lowest level
+ tpart_v(i,kte,j)= tpart_v(i,kte-1,j) ! highest level
+ ENDDO
+ ENDDO
+ ENDIF ! endif (trigger.eq.2)
!
DO J = jts,jte
DO I= its,ite
@@ -282,7 +493,24 @@ SUBROUTINE KF_eta_CPS( &
DQRDT(k)=0.
DQSDT(k)=0.
DTDT(k)=0.
+!ckay
+ cldfra_dp_KF(I,k,J)=0.
+ cldfra_sh_KF(I,k,J)=0.
+ qc_KF(I,k,J)=0.
+ qi_KF(I,k,J)=0.
ENDDO
+#if !(defined(mpas))
+!kf_edrates
+ IF (KF_EDRATES == 1) THEN
+ DO k=kts,kte
+ UDR_KF(I,k,J)=0.
+ DDR_KF(I,k,J)=0.
+ UER_KF(I,k,J)=0.
+ DER_KF(I,k,J)=0.
+ ENDDO
+ TIMEC_KF(I,J)=0.
+ ENDIF
+#endif
RAINCV(I,J)=0.
CUTOP(I,J)=KTS
CUBOT(I,J)=KTE+1
@@ -299,19 +527,47 @@ SUBROUTINE KF_eta_CPS( &
P1D(K) =Pcps(I,K,J)
W0AVG1D(K) =W0AVG(I,K,J)
DZ1D(k)=dz8w(I,K,J)
+ IF (trigger.eq.2) THEN
+ tpart_h1D(K) =tpart_h(I,K,J)
+ tpart_v1D(K) =tpart_v(I,K,J)
+ ELSE
+ tpart_h1D(K) = 0.
+ tpart_v1D(K) = 0.
+ ENDIF
ENDDO
CALL KF_eta_PARA(I, J, &
- U1D,V1D,T1D,QV1D,P1D,DZ1D, &
- W0AVG1D,DT,DX,DXSQ,RHO1D, &
+ U1D,V1D,T1D,QV1D,P1D,DZ1D,W0AVG1D, &
+ tpart_h1D,tpart_v1D, &
+ trigger, &
+ DT,DX,DXSQ,RHO1D, &
XLV0,XLV1,XLS0,XLS1,CP,R,G, &
EP2,SVP1,SVP2,SVP3,SVPT0, &
DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
RAINCV,PRATEC,NCA, &
flag_QI,flag_QS,warm_rain, &
+#if defined(mpas)
CUTOP,CUBOT,CUDT, &
+#else
+ CUTOP,CUBOT,ishall,CUDT, & !BSINGH - Added ishall as arg for WRFCuP
+#endif
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
+ its,ite, jts,jte, kts,kte, &
+!ckay
+ cldfra_dp_KF,cldfra_sh_KF, &
+ qc_KF,qi_KF &
+#if !(defined(mpas))
+!kf_edrates
+ ,UDR_KF,DDR_KF, &
+ UER_KF,DER_KF, &
+ TIMEC_KF,KF_EDRATES &
+#endif
+ )
+#if !(defined(mpas))
+ if(present(shall))shall(i,j) = ishall !,BSINGH - For WRFCuP scheme
+#endif
+
+
IF(PRESENT(rthcuten).AND.PRESENT(rqvcuten)) THEN
DO K=kts,kte
RTHCUTEN(I,K,J)=DTDT(K)/pi(I,K,J)
@@ -356,23 +612,38 @@ SUBROUTINE KF_eta_CPS( &
ENDIF
ENDDO ! i-loop
ENDDO ! j-loop
- ENDIF ! run_param
!
END SUBROUTINE KF_eta_CPS
! ****************************************************************************
!-----------------------------------------------------------
SUBROUTINE KF_eta_PARA (I, J, &
U0,V0,T0,QV0,P0,DZQ,W0AVG1D, &
+ TPART_H0,TPART_V0, &
+ trigger, &
DT,DX,DXSQ,rhoe, &
XLV0,XLV1,XLS0,XLS1,CP,R,G, &
EP2,SVP1,SVP2,SVP3,SVPT0, &
DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
RAINCV,PRATEC,NCA, &
F_QI,F_QS,warm_rain, &
+#if defined(mpas)
CUTOP,CUBOT,CUDT, &
+#else
+ CUTOP,CUBOT,ishall,CUDT, & !BSINGH - For WRFCuP scheme
+#endif
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
+ its,ite, jts,jte, kts,kte, &
+!ckay
+ cldfra_dp_KF,cldfra_sh_KF, &
+ qc_KF,qi_KF &
+#if !(defined(mpas))
+!kf_edrates
+ ,UDR_KF,DDR_KF, &
+ UER_KF,DER_KF, &
+ TIMEC_KF,KF_EDRATES &
+#endif
+ )
!-----------------------------------------------------------
!***** The KF scheme that is currently used in experimental runs of EMCs
!***** Eta model....jsk 8/00
@@ -384,6 +655,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
its,ite, jts,jte, kts,kte, &
I,J
! ,P_QI,P_QS,P_FIRST_SCALAR
+ INTEGER, INTENT(IN ) :: trigger
LOGICAL, INTENT(IN ) :: F_QI, F_QS
@@ -392,6 +664,8 @@ SUBROUTINE KF_eta_PARA (I, J, &
REAL, DIMENSION( kts:kte ), &
INTENT(IN ) :: U0, &
V0, &
+ TPART_H0, &
+ TPART_V0, &
T0, &
QV0, &
P0, &
@@ -406,6 +680,11 @@ SUBROUTINE KF_eta_PARA (I, J, &
REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
!
+#if defined(mpas)
+ INTEGER :: ISHALL ! make ISHALL local for now
+#else
+ INTEGER, INTENT(INOUT) :: ISHALL
+#endif
REAL, DIMENSION( kts:kte ), INTENT(INOUT) :: &
DQDT, &
DQIDT, &
@@ -417,6 +696,26 @@ SUBROUTINE KF_eta_PARA (I, J, &
REAL, DIMENSION( ims:ime , jms:jme ), &
INTENT(INOUT) :: NCA
+!ckay
+ REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
+ INTENT(INOUT) :: cldfra_dp_KF, &
+ cldfra_sh_KF, &
+ qc_KF, &
+ qi_KF
+#if !(defined(mpas))
+!kf_edrates
+ REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
+ INTENT(INOUT) :: UDR_KF, &
+ DDR_KF, &
+ UER_KF, &
+ DER_KF
+
+ REAL, DIMENSION( ims:ime , jms:jme ), &
+ INTENT(INOUT) :: TIMEC_KF
+
+ INTEGER, INTENT(IN) :: KF_EDRATES
+#endif
+
REAL, DIMENSION( ims:ime , jms:jme ), &
INTENT(INOUT) :: RAINCV
@@ -465,14 +764,14 @@ SUBROUTINE KF_eta_PARA (I, J, &
CPORQ,PLCL,ES,DLP,TENV,QENV,TVEN,TVBAR, &
ZLCL,WKL,WABS,TRPPT,WSIGNE,DTLCL,GDT,WLCL,&
TVAVG,QESE,WTW,RHOLCL,AU0,VMFLCL,UPOLD, &
- UPNEW,ABE,WKLCL,TTEMP,FRC1, &
+ UPNEW,ABE,WKLCL,TTEMP,FRC1, &
QNEWIC,RL,R1,QNWFRZ,EFFQ,BE,BOTERM,ENTERM,&
DZZ,UDLBE,REI,EE2,UD2,TTMP,F1,F2, &
THTTMP,QTMP,TMPLIQ,TMPICE,TU95,TU10,EE1, &
UD1,DPTT,QNEWLQ,DUMFDP,EE,TSAT, &
- THTA,VCONV,TIMEC,SHSIGN,VWS,PEF, &
+ THTA,UPFLX,TIMEC,SHSIGN,VWS,PEF, &
CBH,RCBH,PEFCBH,PEFF,PEFF2,TDER,THTMIN, &
- DTMLTD,QS,TADVEC,DPDD,FRC,DPT,RDD,A1, &
+ DTMLTD,QS,TCONV,DPDD,FRC,DPT,RDD,A1, &
DSSDT,DTMP,T1RH,QSRH,PPTFLX,CPR,CNDTNF, &
UPDINC,AINCM2,DEVDMF,PPR,RCED,DPPTDF, &
DMFLFS,DMFLFS2,RCED2,DDINC,AINCMX,AINCM1, &
@@ -488,6 +787,9 @@ SUBROUTINE KF_eta_PARA (I, J, &
INTEGER :: INDLU,NU,NUCHM,NNN,KLFS
REAL :: CHMIN,PM15,CHMAX,DTRH,RAD,DPPP
REAL :: TVDIFF,DTTOT,ABSOMG,ABSOMGTC,FRDP
+ REAL :: DPLUME
+!ckay
+ REAL :: xcldfra,UMF_new
INTEGER :: KX,K,KL
!
@@ -501,16 +803,18 @@ SUBROUTINE KF_eta_PARA (I, J, &
LTOPM1,LVF,KSTART,KMIN,LFS, &
ND,NIC,LDB,LDT,ND1,NDK, &
NM,LMAX,NCOUNT,NOITR, &
- NSTEP,NTC,NCHM,ISHALL,NSHALL
+ NSTEP,NTC,NCHM,NSHALL
LOGICAL :: IPRNT
+ REAL :: u00,qslcl,rhlcl,dqssdt !jfb
CHARACTER*1024 message
!
DATA P00,T00/1.E5,273.16/
DATA RLF/3.339E5/
DATA RHIC,RHBC/1.,0.90/
DATA PIE,TTFRZ,TBFRZ,C5/3.141592654,268.16,248.16,1.0723E-3/
- DATA RATE/0.03/
+ DATA RATE/0.03/ ! wrf default
! DATA RATE/0.01/ ! value used in NRCM
+! DATA RATE/0.001/ ! effectively turn off autoconversion
!-----------------------------------------------------------
IPRNT=.FALSE.
GDRY=-G/CP
@@ -518,6 +822,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
NSHALL = 0
KL=kte
KX=kte
+ TCONV=0.0
!
! ALIQ = 613.3
! BLIQ = 17.502
@@ -556,6 +861,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
!
DO K=1,KX
!
+! Saturation vapor pressure (ES) is calculated following Buck (1981)
!...IF Q0 IS ABOVE SATURATION VALUE, REDUCE IT TO SATURATION LEVEL...
!
ES=ALIQ*EXP((BLIQ*T0(K)-CLIQ)/(T0(K)-DLIQ))
@@ -630,6 +936,13 @@ SUBROUTINE KF_eta_PARA (I, J, &
FBFRC=1.
CYCLE usl
ELSE
+ !BSINGH - For WRFCuP scheme
+ ! wig, 29-Aug-2006: I think this is where no convecion occurs. So, force
+ ! ishall to a flag value to indicate this for accounting purposes.
+#if !(defined(mpas))
+ ishall = 2
+#endif
+ !BSINGH -ENDS
RETURN
ENDIF
ENDIF
@@ -648,13 +961,13 @@ SUBROUTINE KF_eta_PARA (I, J, &
NK=LC-1
IF ( NK+1 .LT. KTS ) THEN
WRITE(message,*)'WOULD GO OFF BOTTOM: KF_ETA_PARA I,J,NK',I,J,NK
-! CALL wrf_message (TRIM(message))
+ WRITE_MESSAGE (TRIM(message))
ELSE
DO
NK=NK+1
IF ( NK .GT. KTE ) THEN
WRITE(message,*)'WOULD GO OFF TOP: KF_ETA_PARA I,J,DPTHMX,DPMIN',I,J,DPTHMX,DPMIN
-! CALL wrf_message (TRIM(message))
+ WRITE_MESSAGE (TRIM(message))
EXIT
ENDIF
DPTHMX=DPTHMX+DP(NK)
@@ -665,6 +978,12 @@ SUBROUTINE KF_eta_PARA (I, J, &
END DO
ENDIF
IF(DPTHMX.LT.DPMIN)THEN
+ !BSINGH - For WRFCuP scheme
+ ! wig, 29-Aug-2006: Indicate no convection occurred in ishall.
+#if !(defined(mpas))
+ ishall = 2
+#endif
+ !BSINGH -ENDS
RETURN
ENDIF
KPBL=LC+NLAYRS-1
@@ -715,18 +1034,36 @@ SUBROUTINE KF_eta_PARA (I, J, &
TLCL=AMIN1(TLCL,TMIX)
TVLCL=TLCL*(1.+0.608*QMIX)
ZLCL = ZMIX+(TLCL-TMIX)/GDRY
- NK = LC-1
- DO
- NK = NK+1
- KLCL=NK
- IF(ZLCL.LE.Z0(NK) .or. NK.GT.KL)THEN
- EXIT
- ENDIF
- ENDDO
- IF(NK.GT.KL)THEN
- RETURN
- ENDIF
+ ! NK = LC-1
+ ! DO
+ ! NK = NK+1
+ ! KLCL=NK
+ ! IF(ZLCL.LE.Z0(NK) .or. NK.GT.KL)THEN
+ ! EXIT
+ ! ENDIF
+ ! ENDDO
+ ! IF(NK.GT.KL)THEN
+ ! RETURN
+ ! ENDIF
+
+ DO NK = LC, KL
+ KLCL = NK
+ IF ( ZLCL.LE.Z0(NK) ) EXIT
+ END DO
+ IF ( ZLCL.GT.Z0(KL) ) then
+ !BSINGH - For WRFCuP scheme
+ !BSINGH - Improvised based on the following reason
+ !Before every "RETURN" statement assign ishall=2
+ ! wig, 29-Aug-2006: Indicate no convection occurred in ishall.
+#if !(defined(mpas))
+ ishall = 2
+#endif
+ !BSINGH -ENDS
+ RETURN
+ endif
+
K=KLCL-1
+! calculate DLP using Z instead of log(P)
DLP=(ZLCL-Z0(K))/(Z0(KLCL)-Z0(K))
!
!...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
@@ -745,41 +1082,51 @@ SUBROUTINE KF_eta_PARA (I, J, &
!...ADJUST VERTICAL VELOCITY TO EQUIVALENT VALUE FOR 25 KM GRID
!...LENGTH, ASSUMING LINEAR DEPENDENCE OF W ON GRID LENGTH...
!
- IF(ZLCL.LT.2.E3)THEN
+ IF(ZLCL.LT.2.E3)THEN ! Kain (2004) Eq. 2
WKLCL=0.02*ZLCL/2.E3
ELSE
- WKLCL=0.02
+ WKLCL=0.02 ! units of m/s
ENDIF
WKL=(W0AVG1D(K)+(W0AVG1D(KLCL)-W0AVG1D(K))*DLP)*DX/25.E3-WKLCL
IF(WKL.LT.0.0001)THEN
DTLCL=0.
ELSE
- DTLCL=4.64*WKL**0.33
+ DTLCL=4.64*WKL**0.33 ! Kain (2004) Eq. 1
ENDIF
+
+! New trigger function
+ IF(trigger.eq.2) then
+ DTLCL=amax1(TPART_H0(KLCL)+TPART_V0(KLCL),0.0)
+ ENDIF
+! end for trigger function
!
+ dtrh = 0.
+ if (trigger .eq. 3) then
!...for ETA model, give parcel an extra temperature perturbation based
!...the threshold RH for condensation (U00)...
+! as described in Narita and Ohmori (2007, 12th Mesoscale Conf.
!
!...for now, just assume U00=0.75...
!...!!!!!! for MM5, SET DTRH = 0. !!!!!!!!
-! U00 = 0.75
-! IF(U00.lt.1.)THEN
-! QSLCL=QES(K)+(QES(KLCL)-QES(K))*DLP
-! RHLCL = QENV/QSLCL
-! DQSSDT = QMIX*(CLIQ-BLIQ*DLIQ)/((TLCL-DLIQ)*(TLCL-DLIQ))
-! IF(RHLCL.ge.0.75 .and. RHLCL.le.0.95)then
-! DTRH = 0.25*(RHLCL-0.75)*QMIX/DQSSDT
-! ELSEIF(RHLCL.GT.0.95)THEN
-! DTRH = (1./RHLCL-1.)*QMIX/DQSSDT
-! ELSE
+ U00 = 0.75
+ IF(U00.lt.1.)THEN
+ QSLCL=QES(K)+(QES(KLCL)-QES(K))*DLP
+ RHLCL = QENV/QSLCL
+ DQSSDT = QMIX*(CLIQ-BLIQ*DLIQ)/((TLCL-DLIQ)*(TLCL-DLIQ))
+ IF(RHLCL.ge.0.75 .and. RHLCL.le.0.95)then
+ DTRH = 0.25*(RHLCL-0.75)*QMIX/DQSSDT
+ ELSEIF(RHLCL.GT.0.95)THEN
+ DTRH = (1./RHLCL-1.)*QMIX/DQSSDT
+ ELSE
DTRH = 0.
-! ENDIF
-! ENDIF
+ ENDIF
+ ENDIF
+ endif ! trigger 3
! IF(ISHALL.EQ.1)IPRNT=.TRUE.
! IPRNT=.TRUE.
! IF(TLCL+DTLCL.GT.TENV)GOTO 45
-!
-trigger: IF(TLCL+DTLCL+DTRH.LT.TENV)THEN
+
+trigger2: IF(TLCL+DTLCL+DTRH.LT.TENV)THEN
!
! Parcel not buoyant, CYCLE back to start of trigger and evaluate next potential USL...
!
@@ -797,7 +1144,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
!
DTTOT = DTLCL+DTRH
IF(DTTOT.GT.1.E-4)THEN
- GDT=2.*G*DTTOT*500./TVEN
+ GDT=2.*G*DTTOT*500./TVEN ! Kain (2004) Eq. 3 (sort of)
WLCL=1.+0.5*SQRT(GDT)
WLCL = AMIN1(WLCL,3.)
ELSE
@@ -811,7 +1158,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
!
LCL=KLCL
LET=LCL
-! make RAD a function of background vertical velocity...
+! make RAD a function of background vertical velocity... (Kain (2004) Eq. 6)
IF(WKL.LT.0.)THEN
RAD = 1000.
ELSEIF(WKL.GT.0.1)THEN
@@ -872,6 +1219,8 @@ SUBROUTINE KF_eta_PARA (I, J, &
!...MIXING RATIO, VERTICAL MASS FLUX, LATERAL DETRAINMENT OF MASS AND
!...MOISTURE, PRECIPITATION RATES AT EACH MODEL LEVEL...
!
+! **1 variables indicate the bottom of a model layer
+! **2 variables indicate the top of a model layer
!
EE1=1.
UD1=0.
@@ -949,7 +1298,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
!
!...REI IS THE RATE OF ENVIRONMENTAL INFLOW...
!
- REI=VMFLCL*DP(NK1)*0.03/RAD
+ REI=VMFLCL*DP(NK1)*0.03/RAD ! KF (1990) Eq. 1; Kain (2004) Eq. 5
TVQU(NK1)=TU(NK1)*(1.+0.608*QU(NK1)-QLIQ(NK1)-QICE(NK1))
IF(NK.EQ.K)THEN
DILBE=((TVLCL+TVQU(NK1))/(TVEN+TV0(NK1))-1.)*DZZ
@@ -962,7 +1311,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
!...ENTRAINMENT (0.5*REI) IS IMPOSED...
!
IF(TVQU(NK1).LE.TV0(NK1))THEN ! Entrain/Detrain IF BLOCK
- EE2=0.5
+ EE2=0.5 ! Kain (2004) Eq. 4
UD2=1.
EQFRC(NK1)=0.
ELSE
@@ -1077,12 +1426,12 @@ SUBROUTINE KF_eta_PARA (I, J, &
!
END DO updraft
!
-!...CHECK CLOUD DEPTH...IF CLOUD IS TALL ENOUGH, ESTIMATE THE EQUILIBRIU
+!...CHECK CLOUD DEPTH...IF CLOUD IS TALL ENOUGH, ESTIMATE THE EQUILIBRIUM
! TEMPERATURE LEVEL (LET) AND ADJUST MASS FLUX PROFILE AT CLOUD TOP SO
-! THAT MASS FLUX DECREASES TO ZERO AS A LINEAR FUNCTION OF PRESSURE BE
+! THAT MASS FLUX DECREASES TO ZERO AS A LINEAR FUNCTION OF PRESSURE BETWEEN
! THE LET AND CLOUD TOP...
!
-!...LTOP IS THE MODEL LEVEL JUST BELOW THE LEVEL AT WHICH VERTICAL VELOC
+!...LTOP IS THE MODEL LEVEL JUST BELOW THE LEVEL AT WHICH VERTICAL VELOCITY
! FIRST BECOMES NEGATIVE...
!
LTOP=NK
@@ -1091,7 +1440,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
!...Instead of using the same minimum cloud height (for deep convection)
!...everywhere, try specifying minimum cloud depth as a function of TLCL...
!
-!
+! Kain (2004) Eq. 7
!
IF(TLCL.GT.293.)THEN
CHMIN = 4.E3
@@ -1100,6 +1449,11 @@ SUBROUTINE KF_eta_PARA (I, J, &
ELSEIF(TLCL.LT.273.)THEN
CHMIN = 2.E3
ENDIF
+!ckay
+ DO NK=K,LTOP
+ qc_KF(I,NK,J)=QLIQ(NK)
+ qi_KF(I,NK,J)=QICE(NK)
+ END DO
!
!...If cloud top height is less than the specified minimum for deep
@@ -1127,15 +1481,28 @@ SUBROUTINE KF_eta_PARA (I, J, &
DETIC(NK)=0.
PPTLIQ(NK)=0.
PPTICE(NK)=0.
+!ckay
+ cldfra_dp_KF(I,NK,J)=0.
+ cldfra_sh_KF(I,NK,J)=0.
+ qc_KF(I,NK,J)=0.
+ qi_KF(I,NK,J)=0.
ENDDO
!
ELSEIF(CLDHGT(LC).GT.CHMIN .and. ABE.GT.1)THEN ! Deep Convection allowed
ISHALL=0
+!ckay
+ DO NK=K,LTOP
+ cldfra_sh_KF(I,NK,J)=0.
+ ENDDO
EXIT usl
ELSE
!
!...TO DISALLOW SHALLOW CONVECTION, COMMENT OUT NEXT LINE !!!!!!!!
ISHALL = 1
+!ckay
+ DO NK=K,LTOP
+ cldfra_dp_KF(I,NK,J)=0.
+ ENDDO
IF(NU.EQ.NUCHM)THEN
EXIT usl ! Shallow Convection from this layer
ELSE
@@ -1148,10 +1515,15 @@ SUBROUTINE KF_eta_PARA (I, J, &
DETIC(NK)=0.
PPTLIQ(NK)=0.
PPTICE(NK)=0.
+!ckay
+ cldfra_dp_KF(I,NK,J)=0.
+ cldfra_sh_KF(I,NK,J)=0.
+ qc_KF(I,NK,J)=0.
+ qi_KF(I,NK,J)=0.
ENDDO
ENDIF
ENDIF
- ENDIF trigger
+ ENDIF trigger2
END DO usl
IF(ISHALL.EQ.1)THEN
KSTART=MAX0(KPBL,KLCL)
@@ -1184,7 +1556,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
!
!...entrainment is allowed at every level except for LTOP, so disallow
!...entrainment at LTOP and adjust entrainment rates between LET and LTOP
-!...so the the dilution factor due to entyrianment is not changed but
+!...so the the dilution factor due to entrainment is not changed but
!...the actual entrainment rate will change due due forced total
!...detrainment in this layer...
!
@@ -1235,6 +1607,11 @@ SUBROUTINE KF_eta_PARA (I, J, &
UMF(NK)=0.
WU(NK)=0.
UER(NK)=0.
+!ckay
+ cldfra_dp_KF(I,NK,J)=0.
+ cldfra_sh_KF(I,NK,J)=0.
+ qc_KF(I,NK,J)=0.
+ qi_KF(I,NK,J)=0.
ENDIF
UDR(NK)=0.
QDT(NK)=0.
@@ -1273,6 +1650,11 @@ SUBROUTINE KF_eta_PARA (I, J, &
TU(NK)=0.
QU(NK)=0.
WU(NK)=0.
+!ckay
+ cldfra_dp_KF(I,NK,J)=0.
+ cldfra_sh_KF(I,NK,J)=0.
+ qc_KF(I,NK,J)=0.
+ qi_KF(I,NK,J)=0.
ENDIF
THTA0(NK)=0.
THTAU(NK)=0.
@@ -1291,7 +1673,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
EMS(NK)=DP(NK)*DXSQ/G
EMSD(NK)=1./EMS(NK)
!
-!...INITIALIZE SOME VARIABLES TO BE USED LATER IN THE VERT ADVECTION SCH
+!...INITIALIZE SOME VARIABLES TO BE USED LATER IN THE VERT ADVECTION SCHEME
!
EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QDT(NK)))
THTAU(NK)=TU(NK)*EXN(NK)
@@ -1307,25 +1689,27 @@ SUBROUTINE KF_eta_PARA (I, J, &
! * WLCL,CLDHGT
! ENDIF
!
-!...COMPUTE CONVECTIVE TIME SCALE(TIMEC). THE MEAN WIND AT THE LCL
-!...AND MIDTROPOSPHERE IS USED.
-!
- WSPD(KLCL)=SQRT(U0(KLCL)*U0(KLCL)+V0(KLCL)*V0(KLCL))
- WSPD(L5)=SQRT(U0(L5)*U0(L5)+V0(L5)*V0(L5))
- WSPD(LTOP)=SQRT(U0(LTOP)*U0(LTOP)+V0(LTOP)*V0(LTOP))
- VCONV=.5*(WSPD(KLCL)+WSPD(L5))
-!...for ETA model, DX is a function of location...
-! TIMEC=DX(I,J)/VCONV
- TIMEC=DX/VCONV
- TADVEC=TIMEC
- TIMEC=AMAX1(1800.,TIMEC)
- TIMEC=AMIN1(3600.,TIMEC)
- IF(ISHALL.EQ.1)TIMEC=2400.
+!...COMPUTE CONVECTIVE TIME SCALE (TIMEC) BASED ON UPDRAFT DYNAMICS
+!
+ IF(ISHALL.NE.1)THEN
+ UPFLX=UMF(KLCL)/DXSQ
+ DPLUME=P0(KLCL)/(R*TU(KLCL))
+ UPFLX=UPFLX/DPLUME
+ TIMEC=CLDHGT(LC)/((UPFLX*ABE)**0.3333)
+ TCONV=TIMEC
+ TIMEC=AMAX1(1800.,TIMEC)
+ TIMEC=AMIN1(18000.,TIMEC)
+ ELSE
+ TIMEC=2400.
+ ENDIF
+
NIC=NINT(TIMEC/DT)
TIMEC=FLOAT(NIC)*DT
!
!...COMPUTE WIND SHEAR AND PRECIPITATION EFFICIENCY.
!
+ WSPD(LTOP)=SQRT(U0(LTOP)*U0(LTOP)+V0(LTOP)*V0(LTOP))
+ WSPD(KLCL)=SQRT(U0(KLCL)*U0(KLCL)+V0(KLCL)*V0(KLCL))
IF(WSPD(LTOP).GT.WSPD(KLCL))THEN
SHSIGN=1.
ELSE
@@ -1358,8 +1742,8 @@ SUBROUTINE KF_eta_PARA (I, J, &
IF(IPRNT)THEN
! WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
WRITE(message,1035)PEF,PEFCBH,LC,LET,WKL,VWS
-! CALL wrf_message( message )
-! call flush(98)
+ WRITE_MESSAGE( message )
+! flush(98)
endif
! WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
!*****************************************************************
@@ -1425,7 +1809,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
RHBAR = RHBAR+RH(ND)*DP(ND)
ENDDO
RHBAR = RHBAR/DPTT
- DMFFRC = 2.*(1.-RHBAR)
+ DMFFRC = 2.*(1.-RHBAR) ! Kain (2004) eq. 11
DPDD = 0.
!...Calculate melting effect
!... first, compute total frozen precipitation generated...
@@ -1542,15 +1926,16 @@ SUBROUTINE KF_eta_PARA (I, J, &
DER(NK)=DER(NK)*DDINC
DDR(NK)=DDR(NK)*DDINC
ENDDO
- CPR=TRPPT
- PPTFLX = TRPPT-TDER
- PEFF=PPTFLX/TRPPT
- IF(IPRNT)THEN
+ CPR=TRPPT
+ PPTFLX = TRPPT-TDER
+ IF(PPTFLX.LT.0.) PPTFLX=0. !JAH
+ PEFF=PPTFLX/TRPPT
+ IF(IPRNT)THEN
! write(98,*)'PRECIP EFFICIENCY =',PEFF
- write(message,*)'PRECIP EFFICIENCY =',PEFF
-! CALL wrf_message(message)
-! call flush(98)
- ENDIF
+ write(message,*)'PRECIP EFFICIENCY =',PEFF
+ WRITE_MESSAGE(message)
+! flush(98)
+ ENDIF
!
!
!...ADJUST UPDRAFT MASS FLUX, MASS DETRAINMENT RATE, AND LIQUID WATER AN
@@ -1597,8 +1982,8 @@ SUBROUTINE KF_eta_PARA (I, J, &
ENDDO
ENDIF d_mf
!
-!...SET LIMITS ON THE UPDRAFT AND DOWNDRAFT MASS FLUXES SO THAT THE INFL
-! INTO CONVECTIVE DRAFTS FROM A GIVEN LAYER IS NO MORE THAN IS AVAILAB
+!...SET LIMITS ON THE UPDRAFT AND DOWNDRAFT MASS FLUXES SO THAT THE INFLOW
+! INTO CONVECTIVE DRAFTS FROM A GIVEN LAYER IS NO MORE THAN IS AVAILABLE
! IN THAT LAYER INITIALLY...
!
AINCMX=1000.
@@ -1709,7 +2094,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
!
DO NTC=1,NSTEP
!
-!...ASSIGN THETA AND Q VALUES AT THE TOP AND BOTTOM OF EACH LAYER BASED
+!...ASSIGN THETA AND Q VALUES AT THE TOP AND BOTTOM OF EACH LAYER BASED ON
!...SIGN OF OMEGA...
!
DO NK=1,LTOP
@@ -1747,7 +2132,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
QG(NK)=QPA(NK)
ENDDO
!
-!...CHECK TO SEE IF MIXING RATIO DIPS BELOW ZERO ANYWHERE; IF SO, BORRO
+!...CHECK TO SEE IF MIXING RATIO DIPS BELOW ZERO ANYWHERE; IF SO, BORROW
!...MOISTURE FROM ADJACENT LAYERS TO BRING IT BACK UP ABOVE ZERO...
!
DO NK=1,LTOP
@@ -1755,7 +2140,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
IF(NK.EQ.1)THEN ! JSK MODS
! PRINT *,' PROBLEM WITH KF SCHEME: ' ! JSK MODS
! PRINT *,'QG = 0 AT THE SURFACE!!!!!!!' ! JSK MODS
-! CALL wrf_error_fatal ( 'QG, QG(NK).LT.0') ! JSK MODS
+ FATAL_ERROR ( 'QG, QG(NK).LT.0') ! JSK MODS
ENDIF ! JSK MODS
NK1=NK+1
IF(NK.EQ.LTOP)THEN
@@ -1913,6 +2298,14 @@ SUBROUTINE KF_eta_PARA (I, J, &
IF(FABE.GT.1. .and. ISHALL.EQ.0)THEN
! WRITE(98,*)'UPDRAFT/DOWNDRAFT COUPLET INCREASES CAPE AT THIS
! *GRID POINT; NO CONVECTION ALLOWED!'
+ !BSINGH - For WRFCuP scheme
+ !BSINGH - Improvised based on the following reason
+ !Before every "RETURN" statement assign ishall=2
+ ! wig, 29-Aug-2006: Indicate no convection occurred in ishall.
+#if !(defined(mpas))
+ ishall = 2
+#endif
+ !BSINGH -ENDS
RETURN
ENDIF
IF(NCOUNT.NE.1)THEN
@@ -1948,7 +2341,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
EXIT
ENDIF
!
-!...IF MORE THAN 10% OF THE ORIGINAL CAPE REMAINS, INCREASE THE CONVECTI
+!...IF MORE THAN 10% OF THE ORIGINAL CAPE REMAINS, INCREASE THE CONVECTIVE
!...MASS FLUX BY THE FACTOR AINC:
!
IF(FABE.EQ.0.)THEN
@@ -1967,13 +2360,22 @@ SUBROUTINE KF_eta_PARA (I, J, &
!...IF AINC BECOMES VERY SMALL, EFFECTS OF CONVECTION ! JSK MODS
!...WILL BE MINIMAL SO JUST IGNORE IT... ! JSK MODS
IF(AINC.LT.0.05)then
+ !BSINGH - For WRFCuP scheme
+ !BSINGH - Improvised based on the following reason
+ !Before every "RETURN" statement assign ishall=2
+ ! wig, 29-Aug-2006: Indicate no convection occurred in ishall.
+#if !(defined(mpas))
+ ishall = 2
+#endif
+ !BSINGH -ENDS
+
RETURN ! JSK MODS
ENDIF
! AINC=AMAX1(AINC,0.05) ! JSK MODS
TDER=TDER2*AINC
PPTFLX=PPTFL2*AINC
! IF (XTIME.LT.10.)THEN
-! WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,
+! WRITE(98,1080)LFS,LDB,LDT,TIMEC,TCONV,NSTEP,NCOUNT,
! * FABEOLD,AINCOLD
! ENDIF
DO NK=1,LTOP
@@ -1991,6 +2393,36 @@ SUBROUTINE KF_eta_PARA (I, J, &
!
ENDIF
ENDDO iter
+!ckay
+! get the cloud fraction for layer NK+1=NK1
+ IF(ISHALL.EQ.1) THEN
+ DO NK=KLCL-1, LTOP1
+ UMF_new = UMF(NK)/DXSQ
+ xcldfra = 0.07*alog(1.+(500.*UMF_new))
+ xcldfra = amax1(0.01,xcldfra)
+ cldfra_sh_KF(I,NK,J) = amin1(0.2,xcldfra)
+ ENDDO
+ ELSE
+ DO NK=KLCL-1, LTOP1
+ UMF_new = UMF(NK)/DXSQ
+ xcldfra = 0.14*alog(1.+(500.*UMF_new))
+ xcldfra = amax1(0.01,xcldfra)
+ cldfra_dp_KF(I,NK,J) = amin1(0.6,xcldfra)
+ ENDDO
+ ENDIF
+
+#if !(defined(mpas))
+!kf_edrates
+!Save up/down entrainment/detrainment rates as 3D variables
+ IF (KF_EDRATES == 1) THEN
+ DO NK=1,LTOP
+ UDR_KF(I,NK,J)=UDR(NK)
+ DDR_KF(I,NK,J)=DDR(NK)
+ UER_KF(I,NK,J)=UER(NK)
+ DER_KF(I,NK,J)=DER(NK)
+ ENDDO
+ ENDIF
+#endif
!
!...COMPUTE HYDROMETEOR TENDENCIES AS IS DONE FOR T, QV...
!
@@ -2014,7 +2446,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
ENDDO
DO NTC=1,NSTEP
!
-!...ASSIGN HYDROMETEORS CONCENTRATIONS AT THE TOP AND BOTTOM OF EACH LAY
+!...ASSIGN HYDROMETEORS CONCENTRATIONS AT THE TOP AND BOTTOM OF EACH LAYER
!...BASED ON THE SIGN OF OMEGA...
!
DO NK=1,LTOP
@@ -2064,18 +2496,26 @@ SUBROUTINE KF_eta_PARA (I, J, &
QRG(NK)=QRPA(NK)
QSG(NK)=QSPA(NK)
ENDDO
+
+#if !(defined(mpas))
+!kf_edrates
+!Save convective timescale (TIMEC) as 2D variable
+ IF (KF_EDRATES == 1) THEN
+ TIMEC_KF(I,J)=TIMEC
+ ENDIF
+#endif
!
!...CLEAN THINGS UP, CALCULATE CONVECTIVE FEEDBACK TENDENCIES FOR THIS
!...GRID POINT...
!
! IF (XTIME.LT.10.)THEN
-! WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
+! WRITE(98,1080)LFS,LDB,LDT,TIMEC,TCONV,NSTEP,NCOUNT,FABE,AINC
! ENDIF
IF(IPRNT)THEN
-! WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
- WRITE(message,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
-! CALL wrf_message(message)
-! call flush(98)
+! WRITE(98,1080)LFS,LDB,LDT,TIMEC,TCONV,NSTEP,NCOUNT,FABE,AINC
+ WRITE(message,1080)LFS,LDB,LDT,TIMEC,TCONV,NSTEP,NCOUNT,FABE,AINC
+ WRITE_MESSAGE(message)
+! flush(98)
endif
!
!...SEND FINAL PARAMETERIZED VALUES TO OUTPUT FILES...
@@ -2088,31 +2528,31 @@ SUBROUTINE KF_eta_PARA (I, J, &
! write(98,*)'At t(h), I, J =',float(NTSD)*72./3600.,I,J
write(message,*)'P(LC), DTP, WKL, WKLCL =',p0(LC)/100., &
TLCL+DTLCL+dtrh-TENV,WKL,WKLCL
-! call wrf_message(message)
+ WRITE_MESSAGE(message)
write(message,*)'TLCL, DTLCL, DTRH, TENV =',TLCL,DTLCL, &
DTRH,TENV
-! call wrf_message(message)
+ WRITE_MESSAGE(message)
WRITE(message,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG, &
TMIX-T00,PMIX,QMIX,ABE
-! call wrf_message(message)
+ WRITE_MESSAGE(message)
WRITE(message,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100., &
WLCL,CLDHGT(LC)
-! call wrf_message(message)
+ WRITE_MESSAGE(message)
WRITE(message,1035)PEF,PEFCBH,LC,LET,WKL,VWS
-! call wrf_message(message)
+ WRITE_MESSAGE(message)
write(message,*)'PRECIP EFFICIENCY =',PEFF
-! call wrf_message(message)
- WRITE(message,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
-! call wrf_message(message)
+ WRITE_MESSAGE(message)
+ WRITE(message,1080)LFS,LDB,LDT,TIMEC,TCONV,NSTEP,NCOUNT,FABE,AINC
+ WRITE_MESSAGE(message)
! ENDIF
!!!!! HERE !!!!!!!
WRITE(message,1070)' P ',' DP ',' DT K/D ',' DR K/D ',' OMG ', &
' DOMGDP ',' UMF ',' UER ',' UDR ',' DMF ',' DER ' &
,' DDR ',' EMS ',' W0 ',' DETLQ ',' DETIC '
-! call wrf_message(message)
+ WRITE_MESSAGE(message)
write(message,*)'just before DO 300...'
-! call wrf_message(message)
-! call flush(98)
+ WRITE_MESSAGE(message)
+! flush(98)
DO NK=1,LTOP
K=LTOP-NK+1
DTT=(TG(K)-T0(K))*86400./TIMEC
@@ -2126,11 +2566,11 @@ SUBROUTINE KF_eta_PARA (I, J, &
UMF(K)/1.E6,UEFRC,UDFRC,DMF(K)/1.E6,DEFRC,DDFRC,EMS(K)/1.E11, &
W0AVG1D(K)*1.E2,DETLQ(K)*TIMEC*EMSD(K)*1.E3,DETIC(K)* &
TIMEC*EMSD(K)*1.E3
-! call wrf_message(message)
+ WRITE_MESSAGE(message)
ENDDO
WRITE(message,1085)'K','P','Z','T0','TG','DT','TU','TD','Q0','QG', &
'DQ','QU','QD','QLG','QIG','QRG','QSG','RH0','RHG'
-! call wrf_message(message)
+ WRITE_MESSAGE(message)
DO NK=1,KL
K=KX-NK+1
DTT=TG(K)-T0(K)
@@ -2150,7 +2590,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
TDC,Q0(K)*1000.,QG(K)*1000.,(QG(K)-Q0(K))*1000.,QU(K)* &
1000.,QD(K)*1000.,QLG(K)*1000.,QIG(K)*1000.,QRG(K)*1000., &
QSG(K)*1000.,RH0,RHG
-! call wrf_message(message)
+ WRITE_MESSAGE(message)
ENDDO
!
!...IF CALCULATIONS ABOVE SHOW AN ERROR IN THE MASS BUDGET, PRINT OUT A
@@ -2173,7 +2613,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
! * U0(K),V0(K),DP(K)/100.,W0AVG(I,J,K)
310 CONTINUE
IF(ISTOP.EQ.1)THEN
-! CALL wrf_error_fatal ( 'KAIN-FRITSCH, istop=1, diags' )
+ FATAL_ERROR ( 'KAIN-FRITSCH, istop=1, diags' )
ENDIF
! ENDIF
4455 format(8f11.3)
@@ -2221,7 +2661,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
WRITE(98,4456)P0(K)/100.,T0(K)-273.16,Q0(K)*1000., &
U0(K),V0(K),W0AVG1D(K),dp(k)/100.,tke(k)
311 CONTINUE
-! call flush(98)
+! flush(98)
! GOTO 297
! STOP 'QVERR'
@@ -2241,10 +2681,10 @@ SUBROUTINE KF_eta_PARA (I, J, &
!
!...FEEDBACK TO RESOLVABLE SCALE TENDENCIES.
!
-!...IF THE ADVECTIVE TIME PERIOD (TADVEC) IS LESS THAN SPECIFIED MINIMUM
-!...TIMEC, ALLOW FEEDBACK TO OCCUR ONLY DURING TADVEC...
+!...IF THE CONVECTIVE TIME PERIOD (TCONV) IS LESS THAN SPECIFIED MINIMUM
+!...TIMEC, ALLOW FEEDBACK TO OCCUR ONLY DURING TCONV...
!
- IF(TADVEC.LT.TIMEC)NIC=NINT(TADVEC/DT)
+ IF(TCONV.LT.TIMEC)NIC=NINT(TCONV/DT)
NCA(I,J) = REAL(NIC)*DT
IF(ISHALL.EQ.1)THEN
TIMEC = 2400.
@@ -2255,7 +2695,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
DO K=1,KX
! IF(IMOIST(INEST).NE.2)THEN
!
-!...IF HYDROMETEORS ARE NOT ALLOWED, THEY MUST BE EVAPORATED OR SUBLIMAT
+!...IF HYDROMETEORS ARE NOT ALLOWED, THEY MUST BE EVAPORATED OR SUBLIMATED
!...AND FED BACK AS VAPOR, ALONG WITH ASSOCIATED CHANGES IN TEMPERATURE.
!...NOTE: THIS WILL INTRODUCE CHANGES IN THE CONVECTIVE TEMPERATURE AND
!...WATER VAPOR FEEDBACK TENDENCIES AND MAY LEAD TO SUPERSATURATED VALUE
@@ -2274,7 +2714,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
!
!...IF ICE PHASE IS NOT ALLOWED, MELT ALL FROZEN HYDROMETEORS...
!
- IF(.NOT. F_QI .and. warm_rain)THEN
+ IF(warm_rain)THEN
CPM=CP*(1.+0.887*QG(K))
TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
@@ -2282,9 +2722,9 @@ SUBROUTINE KF_eta_PARA (I, J, &
DQIDT(K)=0.
DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
DQSDT(K)=0.
- ELSEIF(.NOT. F_QI .and. .not. warm_rain)THEN
+ ELSEIF(.NOT. F_QS)THEN
!
-!...IF ICE PHASE IS ALLOWED, BUT MIXED PHASE IS NOT, MELT FROZEN HYDROME
+!...IF ICE PHASE IS ALLOWED, BUT MIXED PHASE IS NOT, MELT FROZEN HYDROMETEORS
!...BELOW THE MELTING LEVEL, FREEZE LIQUID WATER ABOVE THE MELTING LEVEL
!
CPM=CP*(1.+0.887*QG(K))
@@ -2297,22 +2737,22 @@ SUBROUTINE KF_eta_PARA (I, J, &
DQIDT(K)=0.
DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
DQSDT(K)=0.
- ELSEIF(F_QI) THEN
+ ELSEIF(F_QS) THEN
!
-!...IF MIXED PHASE HYDROMETEORS ARE ALLOWED, FEED BACK CONVECTIVE TENDEN
+!...IF MIXED PHASE HYDROMETEORS ARE ALLOWED, FEED BACK CONVECTIVE TENDENCIES
!...OF HYDROMETEORS DIRECTLY...
!
DQCDT(K)=(QLG(K)-QL0(K))/TIMEC
- DQIDT(K)=(QIG(K)-QI0(K))/TIMEC
+ DQSDT(K)=(QSG(K)-QS0(K))/TIMEC
DQRDT(K)=(QRG(K)-QR0(K))/TIMEC
- IF (F_QS) THEN
- DQSDT(K)=(QSG(K)-QS0(K))/TIMEC
+ IF (F_QI) THEN
+ DQIDT(K)=(QIG(K)-QI0(K))/TIMEC
ELSE
- DQIDT(K)=DQIDT(K)+(QSG(K)-QS0(K))/TIMEC
+ DQSDT(K)=DQSDT(K)+(QIG(K)-QI0(K))/TIMEC
ENDIF
ELSE
! PRINT *,'THIS COMBINATION OF IMOIST, IEXICE, IICE NOT ALLOWED!'
-! CALL wrf_error_fatal ( 'KAIN-FRITSCH, THIS COMBINATION OF IMOIST, IEXICE, IICE NOT ALLOWED' )
+ FATAL_ERROR ( 'KAIN-FRITSCH, THIS MICROPHYSICS CHOICE IS NOT ALLOWED' )
ENDIF
DTDT(K)=(TG(K)-T0(K))/TIMEC
DQDT(K)=(QG(K)-Q0(K))/TIMEC
@@ -2327,7 +2767,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
! write (6,909)I,J,RNC
! WRITE(98,*)'at NTSD =',NTSD,',No. of KF points activated =',
! * NCCNT
-! call flush(98)
+! flush(98)
1000 FORMAT(' ',10A8)
1005 FORMAT(' ',F6.0,2X,F6.4,2X,F7.3,1X,F6.4,2X,4(F6.3,2X),2(F7.3,1X))
1010 FORMAT(' ',' VERTICAL VELOCITY IS NEGATIVE AT ',F4.0,' MB')
@@ -2347,7 +2787,7 @@ SUBROUTINE KF_eta_PARA (I, J, &
! &DEGREE OF STABILIZATION! FABE= ',F6.4)
1070 FORMAT (16A8)
1075 FORMAT (F8.2,3(F8.2),2(F8.3),F8.2,2F8.3,F8.2,6F8.3)
- 1080 FORMAT(2X,'LFS,LDB,LDT =',3I3,' TIMEC, TADVEC, NSTEP=', &
+ 1080 FORMAT(2X,'LFS,LDB,LDT =',3I3,' TIMEC, TCONV, NSTEP=', &
2(1X,F5.0),I3,'NCOUNT, FABE, AINC=',I2,1X,F5.3,F6.2)
1085 FORMAT (A3,16A7,2A8)
1090 FORMAT (I3,F7.2,F7.0,10F7.2,4F7.3,2F8.3)
@@ -2373,7 +2813,7 @@ END SUBROUTINE KF_eta_PARA
SUBROUTINE TPMIX2(p,thes,tu,qu,qliq,qice,qnewlq,qnewic,XLV1,XLV0)
!
! Lookup table variables:
-! INTEGER, PARAMETER :: (KFNT=250,KFNP=220)
+! INTEGER, PARAMETER :: (KFNT=300,KFNP=250)
! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
! REAL, SAVE, DIMENSION(1:200) :: ALU
@@ -2391,7 +2831,7 @@ SUBROUTINE TPMIX2(p,thes,tu,qu,qliq,qice,qnewlq,qnewic,XLV1,XLV0)
!-----------------------------------------------------------------------
!c******** LOOKUP TABLE VARIABLES... ****************************
-! parameter(kfnt=250,kfnp=220)
+! parameter(kfnt=300,kfnp=250)
!c
! COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp),
! * alu(200),rdpr,rdthk,plutop
@@ -2411,13 +2851,25 @@ SUBROUTINE TPMIX2(p,thes,tu,qu,qliq,qice,qnewlq,qnewic,XLV1,XLV0)
!***********************************************************************
!
! scaling the and tt table index
+ IF (IPTB.LT.1) THEN
+ write(98,*)'*** IPTB OUT OF BOUNDS *** IPTB = ',IPTB
+ write(98,*)' setting IPTB = 1'
+ write(98,*)
+ IPTB = 1
+ ENDIF
bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)
tth=(thes-bth)*rdthk
pp =tth-aint(tth)
ithtb=int(tth)+1
- IF(IPTB.GE.220 .OR. IPTB.LE.1 .OR. ITHTB.GE.250 .OR. ITHTB.LE.1)THEN
- write(98,*)'**** OUT OF BOUNDS *********'
-! call flush(98)
+ IF(IPTB.GE.250 .OR. IPTB.LE.1 .OR. ITHTB.GE.300 .OR. ITHTB.LE.1)THEN
+ write(98,*)'**** OUT OF BOUNDS in KF routine lookup table ****'
+ write(98,*)' IPTB = ',iptb,' ITHTB = ',ithtb
+ write(98,*)
+ IF(IPTB.GT.250 .OR. IPTB.LT.1 .OR. ITHTB.GT.300 .OR. ITHTB.LT.1)THEN
+ print *, 'Abort by module_cu_kfeta: Lookup location is outside table array'
+ call exit (1)
+ ENDIF
+! flush(98)
ENDIF
!
t00=ttab(ithtb ,iptb )
@@ -2504,7 +2956,7 @@ SUBROUTINE DTFRZNEW(TU,P,THTEU,QU,QFRZ,QICE,ALIQ,BLIQ,CLIQ,DLIQ)
!...ALLOW THE FREEZING OF LIQUID WATER IN THE UPDRAFT TO PROCEED AS AN
!...APPROXIMATELY LINEAR FUNCTION OF TEMPERATURE IN THE TEMPERATURE RANGE
!...TTFRZ TO TBFRZ...
-!...FOR COLDER TERMPERATURES, FREEZE ALL LIQUID WATER...
+!...FOR COLDER TEMPERATURES, FREEZE ALL LIQUID WATER...
!...THERMODYNAMIC PROPERTIES ARE STILL CALCULATED WITH RESPECT TO LIQUID WATER
!...TO ALLOW THE USE OF LOOKUP TABLE TO EXTRACT TMP FROM THETAE...
!
@@ -2530,7 +2982,12 @@ SUBROUTINE DTFRZNEW(TU,P,THTEU,QU,QFRZ,QICE,ALIQ,BLIQ,CLIQ,DLIQ)
!...TEMPERATURE TO THE SATURATION VALUE...
!
DQEVAP = QS-QU
- QICE = QICE-DQEVAP
+ IF(DQEVAP.GE.QICE)THEN !JAH
+ DQEVAP=QICE !JAH
+ QICE=0. !JAH
+ ELSE !JAH
+ QICE = QICE-DQEVAP
+ ENDIF !JAH
QU = QU+DQEVAP
PII=(1.E5/P)**(0.2854*(1.-0.28*QU))
THTEU=TU*PII*EXP((3374.6525/TU-2.5403)*QU*(1.+0.81*QU))
@@ -2554,13 +3011,7 @@ SUBROUTINE CONDLOAD(QLIQ,QICE,WTW,DZ,BOTERM,ENTERM,RATE,QNEWLQ, &
REAL, INTENT(IN ) :: DZ,BOTERM,ENTERM,RATE
REAL, INTENT(INOUT) :: QLQOUT,QICOUT,WTW,QLIQ,QICE,QNEWLQ,QNEWIC
REAL :: QTOT,QNEW,QEST,G1,WAVG,CONV,RATIO3,OLDQ,RATIO4,DQ,PPTDRG
-
!
-! 9/18/88...THIS PRECIPITATION FALLOUT SCHEME IS BASED ON THE SCHEME US
-! BY OGURA AND CHO (1973). LIQUID WATER FALLOUT FROM A PARCEL IS CAL-
-! CULATED USING THE EQUATION DQ=-RATE*Q*DT, BUT TO SIMULATE A QUASI-
-! CONTINUOUS PROCESS, AND TO ELIMINATE A DEPENDENCY ON VERTICAL
-! RESOLUTION THIS IS EXPRESSED AS Q=Q*EXP(-RATE*DZ).
QTOT=QLIQ+QICE
QNEW=QNEWLQ+QNEWIC
!
@@ -2572,7 +3023,7 @@ SUBROUTINE CONDLOAD(QLIQ,QICE,WTW,DZ,BOTERM,ENTERM,RATE,QNEWLQ, &
G1=WTW+BOTERM-ENTERM-2.*G*DZ*QEST/1.5
IF(G1.LT.0.0)G1=0.
WAVG=0.5*(SQRT(WTW)+SQRT(G1))
- CONV=RATE*DZ/WAVG
+ CONV=RATE*DZ/WAVG ! KF90 Eq. 9
!
! RATIO3 IS THE FRACTION OF LIQUID WATER IN FRESH CONDENSATE, RATIO4 IS
! THE FRACTION OF LIQUID WATER IN THE TOTAL AMOUNT OF CONDENSATE INVOLV
@@ -2584,7 +3035,7 @@ SUBROUTINE CONDLOAD(QLIQ,QICE,WTW,DZ,BOTERM,ENTERM,RATE,QNEWLQ, &
QTOT=QTOT+0.6*QNEW
OLDQ=QTOT
RATIO4=(0.6*QNEWLQ+QLIQ)/(QTOT+1.E-8)
- QTOT=QTOT*EXP(-CONV)
+ QTOT=QTOT*EXP(-CONV) ! KF90 Eq. 9
!
! DETERMINE THE AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE UPDRAFT
! PARCEL AT THIS LEVEL...
@@ -2623,6 +3074,8 @@ SUBROUTINE PROF5(EQ,EE,UD)
! MATHEMATICS SERIES. JUNE, 1964., MAY, 1968.
! JACK KAIN
! 7/6/89
+! Solves for KF90 Eq. 2
+!
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!-----------------------------------------------------------------------
IMPLICIT NONE
@@ -2659,7 +3112,7 @@ END SUBROUTINE PROF5
SUBROUTINE TPMIX2DD(p,thes,ts,qs,i,j)
!
! Lookup table variables:
-! INTEGER, PARAMETER :: (KFNT=250,KFNP=220)
+! INTEGER, PARAMETER :: (KFNT=300,KFNP=250)
! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
! REAL, SAVE, DIMENSION(1:200) :: ALU
@@ -2678,7 +3131,7 @@ SUBROUTINE TPMIX2DD(p,thes,ts,qs,i,j)
!
!******** LOOKUP TABLE VARIABLES (F77 format)... ****************************
-! parameter(kfnt=250,kfnp=220)
+! parameter(kfnt=300,kfnp=250)
!
! COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp), &
! alu(200),rdpr,rdthk,plutop
@@ -2740,6 +3193,7 @@ SUBROUTINE ENVIRTHT(P1,T1,Q1,THT1,ALIQ,BLIQ,CLIQ,DLIQ)
! CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMPERATURE...
!
! NOTE: Calculations for mixed/ice phase no longer used...jsk 8/00
+! For example, KF90 Eq. 10 no longer used
!
EE=Q1*P1/(0.622+Q1)
! TLOG=ALOG(EE/ALIQ)
@@ -2862,7 +3316,7 @@ subroutine kf_lutab(SVP1,SVP2,SVP3,SVPT0)
IMPLICIT NONE
!--------------------------------------------------------------------
! Lookup table variables
-! INTEGER, SAVE, PARAMETER :: KFNT=250,KFNP=220
+! INTEGER, SAVE, PARAMETER :: KFNT=300,KFNP=250
! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
! REAL, SAVE, DIMENSION(1:200) :: ALU
@@ -2884,7 +3338,7 @@ subroutine kf_lutab(SVP1,SVP2,SVP3,SVPT0)
! tolerance for accuracy of temperature
data toler/0.001/
! top pressure (pascals)
- plutop=5000.0
+ plutop=1000.0
! bottom pressure (pascals)
pbot=110000.0