diff --git a/atmos_cubed_sphere/tools/fv_mp_mod.F90 b/atmos_cubed_sphere/tools/fv_mp_mod.F90 index 5cc231d91..032e1aaf4 100644 --- a/atmos_cubed_sphere/tools/fv_mp_mod.F90 +++ b/atmos_cubed_sphere/tools/fv_mp_mod.F90 @@ -137,15 +137,21 @@ module fv_mp_mod END INTERFACE INTERFACE mp_bcst - MODULE PROCEDURE mp_bcst_i4 + MODULE PROCEDURE mp_bcst_i MODULE PROCEDURE mp_bcst_r4 MODULE PROCEDURE mp_bcst_r8 + MODULE PROCEDURE mp_bcst_1d_r4 + MODULE PROCEDURE mp_bcst_1d_r8 + MODULE PROCEDURE mp_bcst_2d_r4 + MODULE PROCEDURE mp_bcst_2d_r8 MODULE PROCEDURE mp_bcst_3d_r4 MODULE PROCEDURE mp_bcst_3d_r8 MODULE PROCEDURE mp_bcst_4d_r4 MODULE PROCEDURE mp_bcst_4d_r8 - MODULE PROCEDURE mp_bcst_3d_i8 - MODULE PROCEDURE mp_bcst_4d_i8 + MODULE PROCEDURE mp_bcst_1d_i + MODULE PROCEDURE mp_bcst_2d_i + MODULE PROCEDURE mp_bcst_3d_i + MODULE PROCEDURE mp_bcst_4d_i END INTERFACE INTERFACE mp_reduce_min @@ -158,7 +164,7 @@ module fv_mp_mod MODULE PROCEDURE mp_reduce_max_r4 MODULE PROCEDURE mp_reduce_max_r8_1d MODULE PROCEDURE mp_reduce_max_r8 - MODULE PROCEDURE mp_reduce_max_i4 + MODULE PROCEDURE mp_reduce_max_i END INTERFACE INTERFACE mp_reduce_sum @@ -2087,14 +2093,14 @@ end subroutine mp_gather_3d_r8 !------------------------------------------------------------------------------- ! vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv ! ! -! mp_bcst_i4 :: Call SPMD broadcast +! mp_bcst_i :: Call SPMD broadcast ! - subroutine mp_bcst_i4(q) + subroutine mp_bcst_i(q) integer, intent(INOUT) :: q call MPI_BCAST(q, 1, MPI_INTEGER, masterproc, commglobal, ierror) - end subroutine mp_bcst_i4 + end subroutine mp_bcst_i ! ! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! !------------------------------------------------------------------------------- @@ -2129,6 +2135,70 @@ end subroutine mp_bcst_r8 ! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! !------------------------------------------------------------------------------- +!------------------------------------------------------------------------------- +! vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv ! +! +! mp_bcst_1d_r4 :: Call SPMD broadcast +! + subroutine mp_bcst_1d_r4(q, idim) + integer, intent(IN) :: idim + real(kind=4), intent(INOUT) :: q(idim) + + call MPI_BCAST(q, idim, MPI_REAL, masterproc, commglobal, ierror) + + end subroutine mp_bcst_1d_r4 +! +! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +! vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv ! +! +! mp_bcst_1d_r8 :: Call SPMD broadcast +! + subroutine mp_bcst_1d_r8(q, idim) + integer, intent(IN) :: idim + real(kind=8), intent(INOUT) :: q(idim) + + call MPI_BCAST(q, idim, MPI_DOUBLE_PRECISION, masterproc, commglobal, ierror) + + end subroutine mp_bcst_1d_r8 +! +! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +! vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv ! +! +! mp_bcst_2d_r4 :: Call SPMD broadcast +! + subroutine mp_bcst_2d_r4(q, idim, jdim) + integer, intent(IN) :: idim, jdim + real(kind=4), intent(INOUT) :: q(idim,jdim) + + call MPI_BCAST(q, idim*jdim, MPI_REAL, masterproc, commglobal, ierror) + + end subroutine mp_bcst_2d_r4 +! +! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +! vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv ! +! +! mp_bcst_2d_r8 :: Call SPMD broadcast +! + subroutine mp_bcst_2d_r8(q, idim, jdim) + integer, intent(IN) :: idim, jdim + real(kind=8), intent(INOUT) :: q(idim,jdim) + + call MPI_BCAST(q, idim*jdim, MPI_DOUBLE_PRECISION, masterproc, commglobal, ierror) + + end subroutine mp_bcst_2d_r8 +! +! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! +!------------------------------------------------------------------------------- + !------------------------------------------------------------------------------- ! vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv ! ! @@ -2196,15 +2266,15 @@ end subroutine mp_bcst_4d_r8 !------------------------------------------------------------------------------- ! vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv ! ! -! mp_bcst_3d_i8 :: Call SPMD broadcast +! mp_bcst_3d_i :: Call SPMD broadcast ! - subroutine mp_bcst_3d_i8(q, idim, jdim, kdim) + subroutine mp_bcst_3d_i(q, idim, jdim, kdim) integer, intent(IN) :: idim, jdim, kdim integer, intent(INOUT) :: q(idim,jdim,kdim) call MPI_BCAST(q, idim*jdim*kdim, MPI_INTEGER, masterproc, commglobal, ierror) - end subroutine mp_bcst_3d_i8 + end subroutine mp_bcst_3d_i ! ! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! !------------------------------------------------------------------------------- @@ -2212,15 +2282,46 @@ end subroutine mp_bcst_3d_i8 !------------------------------------------------------------------------------- ! vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv ! ! -! mp_bcst_4d_i8 :: Call SPMD broadcast +! mp_bcst_1d_i :: Call SPMD broadcast +! + subroutine mp_bcst_1d_i(q, idim) + integer, intent(IN) :: idim + integer, intent(INOUT) :: q(idim) + + call MPI_BCAST(q, idim, MPI_INTEGER, masterproc, commglobal, ierror) + + end subroutine mp_bcst_1d_i +! +! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +! vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv ! +! +! mp_bcst_2d_i :: Call SPMD broadcast +! + subroutine mp_bcst_2d_i(q, idim, jdim) + integer, intent(IN) :: idim, jdim + integer, intent(INOUT) :: q(idim,jdim) + + call MPI_BCAST(q, idim*jdim, MPI_INTEGER, masterproc, commglobal, ierror) + + end subroutine mp_bcst_2d_i +! +! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! +!------------------------------------------------------------------------------- +!------------------------------------------------------------------------------- +! vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv ! +! +! mp_bcst_4d_i :: Call SPMD broadcast ! - subroutine mp_bcst_4d_i8(q, idim, jdim, kdim, ldim) + subroutine mp_bcst_4d_i(q, idim, jdim, kdim, ldim) integer, intent(IN) :: idim, jdim, kdim, ldim integer, intent(INOUT) :: q(idim,jdim,kdim,ldim) call MPI_BCAST(q, idim*jdim*kdim*ldim, MPI_INTEGER, masterproc, commglobal, ierror) - end subroutine mp_bcst_4d_i8 + end subroutine mp_bcst_4d_i ! ! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! !------------------------------------------------------------------------------- @@ -2334,9 +2435,9 @@ end subroutine mp_reduce_min_r8 !------------------------------------------------------------------------------- ! vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv ! ! -! mp_bcst_4d_i4 :: Call SPMD REDUCE_MAX +! mp_bcst_4d_i :: Call SPMD REDUCE_MAX ! - subroutine mp_reduce_max_i4(mymax) + subroutine mp_reduce_max_i(mymax) integer, intent(INOUT) :: mymax integer :: gmax @@ -2346,7 +2447,7 @@ subroutine mp_reduce_max_i4(mymax) mymax = gmax - end subroutine mp_reduce_max_i4 + end subroutine mp_reduce_max_i ! ! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! !------------------------------------------------------------------------------- diff --git a/atmos_cubed_sphere/tools/fv_nggps_diag.F90 b/atmos_cubed_sphere/tools/fv_nggps_diag.F90 index a94f9590e..ca172b20a 100644 --- a/atmos_cubed_sphere/tools/fv_nggps_diag.F90 +++ b/atmos_cubed_sphere/tools/fv_nggps_diag.F90 @@ -362,6 +362,13 @@ subroutine fv_nggps_diag(Atm, zvir, Time) !--- PS ! Re-compute pressure (dry_mass + water_vapor) surface pressure if(id_ps > 0) then + do k=1,npzo + do j=jsco,jeco + do i=isco,ieco + wk(i,j,k) = Atm(n)%delp(i,j,k)*(1.-sum(Atm(n)%q(i,j,k,2:Atm(n)%flagstruct%nwat))) + enddo + enddo + enddo do j=jsco,jeco do i=isco,ieco psurf(i,j) = ptop diff --git a/gfsphysics/GFS_layer/GFS_physics_driver.F90 b/gfsphysics/GFS_layer/GFS_physics_driver.F90 index 535970fd8..843f691f0 100644 --- a/gfsphysics/GFS_layer/GFS_physics_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_physics_driver.F90 @@ -2897,15 +2897,11 @@ subroutine GFS_physics_driver & Diag%rain(:) = Diag%rainc(:) + frain * rain1(:) if (Model%cal_pre) then ! hchuang: add dominant precipitation type algorithm - i = min(3,Model%num_p3d) ! -! rsun: when ncld = 2 NEED ATTENTION HERE need to re-write this routine -! - call calpreciptype (kdt, Model%nrcm, im, ix, levs, levs+1, & - Tbd%rann, Grid%xlat, Grid%xlon, Stateout%gt0, & - Stateout%gq0, Statein%prsl, Statein%prsi, & - Diag%rain, Statein%phii, Model%num_p3d, & - Sfcprop%tsfc, Diag%sr, Tbd%phy_f3d(1,1,i), & ! input + call calpreciptype (kdt, Model%nrcm, im, ix, levs, levs+1, & + Tbd%rann, Grid%xlat, Grid%xlon, Stateout%gt0, & + Stateout%gq0, Statein%prsl, Statein%prsi, & + Diag%rain, Statein%phii, Sfcprop%tsfc, & !input domr, domzr, domip, doms) ! output ! ! if (lprnt) print*,'debug calpreciptype: DOMR,DOMZR,DOMIP,DOMS ' @@ -2917,14 +2913,18 @@ subroutine GFS_physics_driver & ! & DOMR(i),DOMZR(i),DOMIP(i),DOMS(i) ! end do ! HCHUANG: use new precipitation type to decide snow flag for LSM snow accumulation + + if (Model%imp_physics /= 11) then + do i=1,im + Sfcprop%tprcp(i) = max(0.0, Diag%rain(i) ) + if(doms(i) > 0.0 .or. domip(i) > 0.0) then + Sfcprop%srflag(i) = 1. + else + Sfcprop%srflag(i) = 0. + end if + enddo + endif - do i=1,im - if(doms(i) > 0.0 .or. domip(i) > 0.0) then - Sfcprop%srflag(i) = 1. - else - Sfcprop%srflag(i) = 0. - end if - enddo endif if (Model%lssav) then @@ -2932,6 +2932,13 @@ subroutine GFS_physics_driver & Diag%totice (:) = Diag%totice (:) + Diag%ice(:) Diag%totsnw (:) = Diag%totsnw (:) + Diag%snow(:) Diag%totgrp (:) = Diag%totgrp (:) + Diag%graupel(:) +! + if (Model%cal_pre) then + Diag%tdomr(:) = Diag%tdomr(:) + domr(:) * dtf + Diag%tdomzr(:) = Diag%tdomzr(:) + domzr(:) * dtf + Diag%tdomip(:) = Diag%tdomip(:) + domip(:) * dtf + Diag%tdoms(:) = Diag%tdoms(:) + doms(:) * dtf + endif if (Model%ldiag3d) then Diag%dt3dt(:,:,6) = Diag%dt3dt(:,:,6) + (Stateout%gt0(:,:)-dtdt(:,:)) * frain @@ -2953,51 +2960,32 @@ subroutine GFS_physics_driver & enddo enddo -! --- ... lu: snow-rain detection is performed in land/sice module - - if (Model%cal_pre) then ! hchuang: new precip type algorithm defines srflag - Sfcprop%tprcp(:) = max(0.0, Diag%rain(:)) ! clu: rain -> tprcp -!rsun if (Model%lgfdlmp) then -!rsun do i = 1, im -!rsun Sfcprop%srflag(i) = 0. ! clu: default srflag as 'rain' (i.e. 0) -!rsun! determine convective rain/snow by surface temperature -!rsun! determine large-scale rain/snow by rain/snow coming out directly from MP -!rsun if (Sfcprop%tsfc(i) .ge. 273.15) then -!rsun crain = Diag%rainc(i) -!rsun csnow = 0.0 -!rsun else -!rsun crain = 0.0 -!rsun csnow = Diag%rainc(i) -!rsun endif -!rsun if ((snow0(i,1)+ice0(i,1)+graupel0(i,1)+csnow) .gt. (rain0(i,1)+crain)) then -!rsun Sfcprop%srflag(i) = 1. ! clu: set srflag to 'snow' (i.e. 1) -!rsun endif -!rsun enddo -!rsun endif - else + if (Model%imp_physics == 11) then +! determine convective rain/snow by surface temperature +! determine large-scale rain/snow by rain/snow coming out directly from MP do i = 1, im Sfcprop%tprcp(i) = max(0.0, Diag%rain(i) )! clu: rain -> tprcp Sfcprop%srflag(i) = 0. ! clu: default srflag as 'rain' (i.e. 0) - if (Model%imp_physics == 11) then -! determine convective rain/snow by surface temperature -! determine large-scale rain/snow by rain/snow coming out directly from MP - if (Sfcprop%tsfc(i) .ge. 273.15) then - crain = Diag%rainc(i) - csnow = 0.0 - else - crain = 0.0 - csnow = Diag%rainc(i) - endif - if ((snow0(i,1)+ice0(i,1)+graupel0(i,1)+csnow) .gt. (rain0(i,1)+crain)) then - Sfcprop%srflag(i) = 1. ! clu: set srflag to 'snow' (i.e. 1) - endif + if (Sfcprop%tsfc(i) .ge. 273.15) then + crain = Diag%rainc(i) + csnow = 0.0 else - if (t850(i) <= 273.16) then - Sfcprop%srflag(i) = 1. ! clu: set srflag to 'snow' (i.e. 1) - endif + crain = 0.0 + csnow = Diag%rainc(i) endif - enddo - endif + if ((snow0(i,1)+ice0(i,1)+graupel0(i,1)+csnow) .gt. (rain0(i,1)+crain)) then + Sfcprop%srflag(i) = 1. ! clu: set srflag to 'snow' (i.e. 1) + endif + enddo + else if( .not. Model%cal_pre) then + do i = 1, im + Sfcprop%tprcp(i) = max(0.0, Diag%rain(i) )! clu: rain -> tprcp + Sfcprop%srflag(i) = 0. ! clu: default srflag as 'rain' (i.e. 0) + if (t850(i) <= 273.16) then + Sfcprop%srflag(i) = 1. ! clu: set srflag to 'snow' (i.e. 1) + endif + enddo + endif ! --- ... coupling insertion diff --git a/gfsphysics/GFS_layer/GFS_typedefs.F90 b/gfsphysics/GFS_layer/GFS_typedefs.F90 index 588883b1f..b9916aa33 100644 --- a/gfsphysics/GFS_layer/GFS_typedefs.F90 +++ b/gfsphysics/GFS_layer/GFS_typedefs.F90 @@ -842,6 +842,10 @@ module GFS_typedefs real (kind=kind_phys), pointer :: smcref2(:) => null() !< soil moisture threshold (volumetric) real (kind=kind_phys), pointer :: wet1 (:) => null() !< normalized soil wetness real (kind=kind_phys), pointer :: sr (:) => null() !< snow ratio : ratio of snow to total precipitation + real (kind=kind_phys), pointer :: tdomr (:) => null() !< accumulated rain type + real (kind=kind_phys), pointer :: tdomzr (:) => null() !< accumulated freezing rain type + real (kind=kind_phys), pointer :: tdomip (:) => null() !< accumulated sleet type + real (kind=kind_phys), pointer :: tdoms (:) => null() !< accumulated snow type real (kind=kind_phys), pointer :: skebu_wts(:,:) => null() !< 10 meater u/v wind speed real (kind=kind_phys), pointer :: skebv_wts(:,:) => null() !< 10 meater u/v wind speed @@ -2628,6 +2632,10 @@ subroutine diag_create (Diag, IM, Model) allocate (Diag%smcref2 (IM)) allocate (Diag%wet1 (IM)) allocate (Diag%sr (IM)) + allocate (Diag%tdomr (IM)) + allocate (Diag%tdomzr (IM)) + allocate (Diag%tdomip (IM)) + allocate (Diag%tdoms (IM)) allocate (Diag%skebu_wts(IM,Model%levs)) allocate (Diag%skebv_wts(IM,Model%levs)) allocate (Diag%sppt_wts(IM,Model%levs)) @@ -2701,7 +2709,6 @@ subroutine diag_phys_zero (Diag, Model, linit) Diag%dvsfc = zero Diag%dtsfc = zero Diag%dqsfc = zero -! Diag%totprcp = zero Diag%gflux = zero Diag%dlwsfc = zero Diag%ulwsfc = zero @@ -2712,7 +2719,6 @@ subroutine diag_phys_zero (Diag, Model, linit) Diag%dugwd = zero Diag%dvgwd = zero Diag%psmean = zero -! Diag%cnvprcp = zero Diag%spfhmin = huge Diag%spfhmax = zero Diag%u10mmax = zero @@ -2723,9 +2729,6 @@ subroutine diag_phys_zero (Diag, Model, linit) Diag%ice = zero Diag%snow = zero Diag%graupel = zero -! Diag%totice = zero -! Diag%totsnw = zero -! Diag%totgrp = zero !--- Out Diag%u10m = zero @@ -2755,6 +2758,10 @@ subroutine diag_phys_zero (Diag, Model, linit) Diag%smcref2 = zero Diag%wet1 = zero Diag%sr = zero + Diag%tdomr = zero + Diag%tdomzr = zero + Diag%tdomip = zero + Diag%tdoms = zero Diag%skebu_wts = zero Diag%skebv_wts = zero Diag%sppt_wts = zero diff --git a/gfsphysics/physics/calpreciptype.f90 b/gfsphysics/physics/calpreciptype.f90 index 7260b974d..5376e44cb 100644 --- a/gfsphysics/physics/calpreciptype.f90 +++ b/gfsphysics/physics/calpreciptype.f90 @@ -1,7 +1,7 @@ subroutine calpreciptype(kdt,nrcm,im,ix,lm,lp1,randomno, & xlat,xlon, & gt0,gq0,prsl,prsi,prec, & !input - phii,n3dfercld,tskin,sr,phy_f3d, & !input + phii,tskin, & !input domr,domzr,domip,doms) !output !$$$ subprogram documentation block @@ -25,16 +25,16 @@ subroutine calpreciptype(kdt,nrcm,im,ix,lm,lp1,randomno, & ! ! declare variables. ! - integer,intent(in) :: kdt,nrcm,im,ix,lm,lp1,n3dfercld + integer,intent(in) :: kdt,nrcm,im,ix,lm,lp1 real,intent(in) :: xlat(im),xlon(im) real,intent(in) :: randomno(ix,nrcm) - real(kind=kind_phys),dimension(im), intent(in) :: prec,sr,tskin - real(kind=kind_phys),dimension(ix,lm), intent(in) :: gt0,gq0,prsl,phy_f3d + real(kind=kind_phys),dimension(im), intent(in) :: prec,tskin + real(kind=kind_phys),dimension(ix,lm), intent(in) :: gt0,gq0,prsl real(kind=kind_phys),dimension(ix,lp1),intent(in) :: prsi,phii real(kind=kind_phys),dimension(im), intent(out) :: domr,domzr,domip,doms integer, dimension(nalg) :: sleet,rain,freezr,snow - real(kind=kind_phys),dimension(lm) :: t,q,pmid,f_rimef + real(kind=kind_phys),dimension(lm) :: t,q,pmid real(kind=kind_phys),dimension(lp1) :: pint,zint real(kind=kind_phys), allocatable :: twet(:),rh(:),td(:) ! @@ -69,7 +69,6 @@ subroutine calpreciptype(kdt,nrcm,im,ix,lm,lp1,randomno, & t(k1) = gt0(i,k) q(k1) = gq0(i,k) pmid(k1) = prsl(i,k) ! pressure in pascals - f_rimef(k1) = phy_f3d(i,k) ! ! compute wet bulb temperature ! @@ -105,9 +104,9 @@ subroutine calpreciptype(kdt,nrcm,im,ix,lm,lp1,randomno, & ! debug print statement ! if (abs(xlon(i)*57.29578-114.0) .lt. 0.2 .and. & ! abs(xlat(i)*57.29578-40.0) .lt. 0.2)then -! print*,'debug in calpreciptype: i,im,lm,lp1,xlon,xlat,prec,tskin,sr,nrcm,randomno,n3dfercld ', & -! i,im,lm,lp1,xlon(i)*57.29578,xlat(i)*57.29578,prec(i),tskin(i),sr(i), & -! nrcm,randomno(i,1:nrcm),n3dfercld +! print*,'debug in calpreciptype: i,im,lm,lp1,xlon,xlat,prec,tskin,nrcm,randomno', & +! i,im,lm,lp1,xlon(i)*57.29578,xlat(i)*57.29578,prec(i),tskin(i),, & +! nrcm,randomno(i,1:nrcm) ! do l=1,lm ! print*,'debug in calpreciptype: l,t,q,p,pint,z,twet', & ! l,t(l),q(l), & @@ -182,18 +181,10 @@ subroutine calpreciptype(kdt,nrcm,im,ix,lm,lp1,randomno, & ! ! explicit algorithm (under 18 not admitted without parent or guardian) - if(n3dfercld == 3) then ! ferrier's scheme - call calwxt_explicit(lm,tskin(i),sr(i),f_rimef,iwx) - snow(5) = mod(iwx,2) - sleet(5) = mod(iwx,4)/2 - freezr(5) = mod(iwx,8)/4 - rain(5) = iwx/8 - else - snow(5) = 0 - sleet(5) = 0 - freezr(5) = 0 - rain(5) = 0 - endif + snow(5) = 0 + sleet(5) = 0 + freezr(5) = 0 + rain(5) = 0 ! call calwxt_dominant(nalg,rain(1),freezr(1),sleet(1), & snow(1),domr(i),domzr(i),domip(i),doms(i)) @@ -1280,66 +1271,6 @@ subroutine calwxt_revised(lm,lp1,t,q,pmid,pint, & return end ! -! - subroutine calwxt_explicit(lm,tskin,sr,f_rimef,iwx) -! -! file: calwxt.f -! written: 24 august 2005, g manikin and b ferrier -! -! routine to compute precipitation type using explicit fields -! from the model microphysics - -!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none -! -! list of variables needed -! parameters: -! -! input: - integer, intent(in) :: lm - real,intent(in) :: tskin, sr - real,intent(in) :: f_rimef(lm) - integer,intent(out) :: iwx - real snow -! real psfc -! -! allocate local storage -! - iwx = 0 - -!gsm the rsm is currently incompatible with this routine -!gsm according to b ferrier, there may be a way to write -!gsm a version of this algorithm to work with the rsm -!gsm microphysics, but it doesn't exist at this time - -! a snow ratio less than 0.5 eliminates snow and sleet -! use the skin temperature to distinguish rain from freezing rain -! note that 2-m temperature may be a better choice if the model -! has a cold bias for skin temperature -! - if (sr < 0.5) then -! surface (skin) potential temperature and temperature. -! psfc=pmid(lm) -! tskin=ths*(psfc/p1000)**capa - - if (tskin < 273.15) then ! freezing rain = 4 - iwx = iwx + 4 - else ! rain = 8 - iwx = iwx + 8 - endif - else -! -! distinguish snow from sleet with the rime factor -! - if(f_rimef(lm) >= 10) then ! sleet = 2 - iwx = iwx + 2 - else - snow = 1 - iwx = iwx + 1 - endif - endif - end -! ! subroutine calwxt_dominant(nalg,rain,freezr,sleet,snow, & & domr,domzr,domip,doms) diff --git a/gfsphysics/physics/module_mp_thompson_gfs.F90 b/gfsphysics/physics/module_mp_thompson_gfs.F90 index f677061b6..8c06d8935 100644 --- a/gfsphysics/physics/module_mp_thompson_gfs.F90 +++ b/gfsphysics/physics/module_mp_thompson_gfs.F90 @@ -4076,24 +4076,21 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & INTEGER:: inu_c real, dimension(15), parameter:: g_ratio = (/24,60,120,210,336, & & 504,720,990,1320,1716,2184,2730,3360,4080,4896/) - real :: Nt_c has_qc = .false. has_qi = .false. has_qs = .false. - if(islmski == 1) then - Nt_c = Nt_cl - else - Nt_c = Nt_co - endif - ! print*,'cal_eff:',islmski,Nt_c, Nt_cl,Nt_co do k = kts, kte rho(k) = 0.622*p1d(k)/(R*t1d(k)*(qv1d(k)+0.622)) rc(k) = MAX(R1, qc1d(k)*rho(k)) - nc(k) = Nt_c + if(islmski == 1) then + nc(k) = Nt_cl + else + nc(k) = Nt_co + endif if (rc(k).gt.R1 .and. nc(k).gt.R2) has_qc = .true. ri(k) = MAX(R1, qi1d(k)*rho(k)) ni(k) = MAX(R2, ni1d(k)*rho(k)) diff --git a/io/FV3GFS_io.F90 b/io/FV3GFS_io.F90 index c75d7a6ed..ee6be2e32 100644 --- a/io/FV3GFS_io.F90 +++ b/io/FV3GFS_io.F90 @@ -2373,7 +2373,7 @@ subroutine gfdl_diag_register(Time, Sfcprop, Cldprop, Gfs_diag, Gfs_grid, Atm_bl idx = idx + 1 Diag(idx)%axes = 2 - Diag(idx)%name = 'totice' + Diag(idx)%name = 'totice_ave' Diag(idx)%desc = 'surface ice precipitation rate' Diag(idx)%unit = 'kg/m**2/s' Diag(idx)%mod_name = 'gfs_phys' @@ -2387,7 +2387,7 @@ subroutine gfdl_diag_register(Time, Sfcprop, Cldprop, Gfs_diag, Gfs_grid, Atm_bl idx = idx + 1 Diag(idx)%axes = 2 - Diag(idx)%name = 'totsnw' + Diag(idx)%name = 'totsnw_ave' Diag(idx)%desc = 'surface snow precipitation rate' Diag(idx)%unit = 'kg/m**2/s' Diag(idx)%mod_name = 'gfs_phys' @@ -2401,7 +2401,7 @@ subroutine gfdl_diag_register(Time, Sfcprop, Cldprop, Gfs_diag, Gfs_grid, Atm_bl idx = idx + 1 Diag(idx)%axes = 2 - Diag(idx)%name = 'totgrp' + Diag(idx)%name = 'totgrp_ave' Diag(idx)%desc = 'surface graupel precipitation rate' Diag(idx)%unit = 'kg/m**2/s' Diag(idx)%mod_name = 'gfs_phys' @@ -2732,6 +2732,62 @@ subroutine gfdl_diag_register(Time, Sfcprop, Cldprop, Gfs_diag, Gfs_grid, Atm_bl Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%sr(:) enddo + idx = idx + 1 + Diag(idx)%axes = 2 + Diag(idx)%name = 'crain_ave' + Diag(idx)%desc = 'averaged categorical rain' + Diag(idx)%unit = 'number' + Diag(idx)%mod_name = 'gfs_phys' + Diag(idx)%intpl_method = 'bilinear' + Diag(idx)%cnvfac = cn_one + Diag(idx)%time_avg = .TRUE. + allocate (Diag(idx)%data(nblks)) + do nb = 1,nblks + Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%tdomr(:) + enddo + + idx = idx + 1 + Diag(idx)%axes = 2 + Diag(idx)%name = 'csnow_ave' + Diag(idx)%desc = 'averaged categorical snow' + Diag(idx)%unit = 'number' + Diag(idx)%mod_name = 'gfs_phys' + Diag(idx)%intpl_method = 'bilinear' + Diag(idx)%cnvfac = cn_one + Diag(idx)%time_avg = .TRUE. + allocate (Diag(idx)%data(nblks)) + do nb = 1,nblks + Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%tdoms(:) + enddo + + idx = idx + 1 + Diag(idx)%axes = 2 + Diag(idx)%name = 'cfrzr_ave' + Diag(idx)%desc = 'averaged categorical freezing rain' + Diag(idx)%unit = 'number' + Diag(idx)%mod_name = 'gfs_phys' + Diag(idx)%intpl_method = 'bilinear' + Diag(idx)%cnvfac = cn_one + Diag(idx)%time_avg = .TRUE. + allocate (Diag(idx)%data(nblks)) + do nb = 1,nblks + Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%tdomzr(:) + enddo + + idx = idx + 1 + Diag(idx)%axes = 2 + Diag(idx)%name = 'cicep_ave' + Diag(idx)%desc = 'averaged categorical sleet' + Diag(idx)%unit = 'number' + Diag(idx)%mod_name = 'gfs_phys' + Diag(idx)%intpl_method = 'bilinear' + Diag(idx)%cnvfac = cn_one + Diag(idx)%time_avg = .TRUE. + allocate (Diag(idx)%data(nblks)) + do nb = 1,nblks + Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%tdomip(:) + enddo + idx = idx + 1 Diag(idx)%axes = 3 Diag(idx)%name = 'skebu_wts' @@ -3040,7 +3096,7 @@ subroutine gfdl_diag_register(Time, Sfcprop, Cldprop, Gfs_diag, Gfs_grid, Atm_bl idx = idx + 1 Diag(idx)%axes = 2 Diag(idx)%name = 'crain' - Diag(idx)%desc = 'categorical rain' + Diag(idx)%desc = 'instantaneous categorical rain' Diag(idx)%unit = 'number' Diag(idx)%mod_name = 'gfs_sfc' Diag(idx)%cnvfac = cn_one @@ -3983,7 +4039,6 @@ subroutine fv_phys_bundle_setup(axes, phys_bundle, fcst_grid, quilting, nbdlphys line=__LINE__, & file=__FILE__)) & return ! bail out - if(mpp_pe()==mpp_root_pe())print *,'in fv_phys bundle,i=',i,'physbdl_name=',physbdl_name idx = index(physbdl_name,'_bilinear') if(idx > 0) then outputfile(ibdl) = physbdl_name(1:idx-1) diff --git a/stochastic_physics/get_stochy_pattern.F90 b/stochastic_physics/get_stochy_pattern.F90 index 238c30837..b3534540e 100644 --- a/stochastic_physics/get_stochy_pattern.F90 +++ b/stochastic_physics/get_stochy_pattern.F90 @@ -6,7 +6,7 @@ module get_stochy_pattern_mod use stochy_gg_def use stochy_patterngenerator use stochy_internal_state_mod - use fv_mp_mod, only : mp_reduce_sum,is_master,mp_bcst + use fv_mp_mod, only : mp_reduce_sum,is_master use fv_arrays_mod, only: fv_atmos_type use GFS_typedefs, only: GFS_control_type, GFS_grid_type use mersenne_twister, only: random_seed @@ -14,7 +14,7 @@ module get_stochy_pattern_mod private public get_random_pattern_fv3,get_random_pattern_fv3_vect - public dump_patterns,restore_patterns + public dump_patterns logical :: first_call=.true. #include "mpif.h" contains @@ -305,7 +305,7 @@ end subroutine scalarspect_to_gaugrid subroutine dump_patterns(sfile) implicit none - character*9 :: sfile + character*120 :: sfile integer :: stochlun,k,n stochlun=99 if (is_master()) then @@ -333,70 +333,6 @@ subroutine dump_patterns(sfile) endif close(stochlun) end subroutine dump_patterns - subroutine restore_patterns(sfile) - implicit none - character*9 :: sfile - integer :: stochlun,k,n - stochlun=99 - if (is_master()) then - if (nsppt > 0 .OR. nshum > 0 .OR. nskeb > 0) then - OPEN(stochlun,file=sfile,form='unformatted',status='old') - endif - endif - if (nsppt > 0) then - do n=1,nsppt - call read_pattern(rpattern_sppt(n),1,stochlun) - enddo - endif - if (nshum > 0) then - do n=1,nshum - call read_pattern(rpattern_shum(n),1,stochlun) - enddo - endif - if (nskeb > 0) then - do n=1,nskeb - do k=1,skeblevs - call read_pattern(rpattern_skeb(n),k,stochlun) - enddo - enddo - endif - close(stochlun) - end subroutine restore_patterns -subroutine read_pattern(rpattern,k,lunptn) - implicit none - type(random_pattern), intent(inout) :: rpattern - integer, intent(in) :: lunptn,k - real(kind_dbl_prec),allocatable :: pattern2d(:) - integer nm,nn,ierr - - allocate(pattern2d(2*ndimspec)) - - ! read only on root process, and send to all tasks - if (is_master()) then - read(lunptn) pattern2d - print*,'reading in random pattern (min/max/size)',minval(pattern2d),maxval(pattern2d),size(pattern2d) - endif - do nn=1,2*ndimspec - call mp_bcst(pattern2d(nn)) - enddo - ! subset - do nn=1,len_trie_ls - nm = rpattern%idx_e(nn) - if (nm == 0) cycle - rpattern%spec_e(nn,1,k) = pattern2d(nm) - rpattern%spec_e(nn,2,k) = pattern2d(ndimspec+nm) - enddo - do nn=1,len_trio_ls - nm = rpattern%idx_o(nn) - if (nm == 0) cycle - rpattern%spec_o(nn,1,k) = pattern2d(nm) - rpattern%spec_o(nn,2,k) = pattern2d(ndimspec+nm) - enddo - !print*,'after scatter...',me,maxval(rpattern%spec_e),maxval(rpattern%spec_o) & - ! ,minval(rpattern%spec_e),minval(rpattern%spec_o) - deallocate(pattern2d) - end subroutine read_pattern - subroutine write_pattern(rpattern,lev,lunptn) implicit none type(random_pattern), intent(inout) :: rpattern diff --git a/stochastic_physics/stochastic_physics.F90 b/stochastic_physics/stochastic_physics.F90 index b2d6f3d4c..ae9949fbb 100644 --- a/stochastic_physics/stochastic_physics.F90 +++ b/stochastic_physics/stochastic_physics.F90 @@ -170,7 +170,7 @@ subroutine run_stochastic_physics(nblks,Model,Grid,Coupling) if (do_sppt.EQ. .false. .AND. do_shum.EQ. .false. .and. do_skeb.EQ. .false.) return ! check to see if it is time to write out random patterns if (Model%phour .EQ. fhstoch) then - write(STRFH,FMT='(I6)') nint(Model%fhour) + write(STRFH,FMT='(I6.6)') nint(Model%fhour) sfile='stoch_out.F'//trim(STRFH) call dump_patterns(sfile) endif diff --git a/stochastic_physics/stochy_data_mod.F90 b/stochastic_physics/stochy_data_mod.F90 index 2f5ec1067..d4ec46820 100644 --- a/stochastic_physics/stochy_data_mod.F90 +++ b/stochastic_physics/stochy_data_mod.F90 @@ -107,10 +107,10 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit) ! no spinup needed if initial patterns are defined correctly. spinup_efolds = 0 if (nsppt > 0) then + if (is_master()) print *, 'Initialize random pattern for SPPT' + call patterngenerator_init(sppt_lscale,delt,sppt_tau,sppt,iseed_sppt,rpattern_sppt, & + lonf,latg,jcap,gis_stochy%ls_node,nsppt,1,0) do n=1,nsppt - if (is_master()) print *, 'Initialize random pattern for SPPT',n - call patterngenerator_init(sppt_lscale,delt,sppt_tau,sppt,iseed_sppt,rpattern_sppt, & - lonf,latg,jcap,gis_stochy%ls_node,nsppt,1,0) nspinup = spinup_efolds*sppt_tau(n)/delt if (stochini) then call read_pattern(rpattern_sppt(n),1,stochlun) @@ -140,9 +140,9 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit) endif if (nshum > 0) then if (is_master()) print *, 'Initialize random pattern for SHUM' + call patterngenerator_init(shum_lscale,delt,shum_tau,shum,iseed_shum,rpattern_shum, & + lonf,latg,jcap,gis_stochy%ls_node,nshum,1,0) do n=1,nshum - call patterngenerator_init(shum_lscale,delt,shum_tau,shum,iseed_shum,rpattern_shum, & - lonf,latg,jcap,gis_stochy%ls_node,nshum,1,0) nspinup = spinup_efolds*shum_tau(n)/delt if (stochini) then call read_pattern(rpattern_shum(n),1,stochlun) @@ -176,9 +176,9 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit) skeblevs=nint(skeb_tau(1)/delt*skeb_vdof) ! backscatter noise. if (is_master()) print *, 'Initialize random pattern for SKEB',skeblevs + call patterngenerator_init(skeb_lscale,delt,skeb_tau,skeb,iseed_skeb,rpattern_skeb, & + lonf,latg,jcap,gis_stochy%ls_node,nskeb,skeblevs,skeb_varspect_opt) do n=1,nskeb - call patterngenerator_init(skeb_lscale,delt,skeb_tau,skeb,iseed_skeb,rpattern_skeb, & - lonf,latg,jcap,gis_stochy%ls_node,nskeb,skeblevs,skeb_varspect_opt) do k=1,skeblevs nspinup = spinup_efolds*skeb_tau(n)/delt if (stochini) then @@ -303,12 +303,8 @@ subroutine read_pattern(rpattern,k,lunptn) endif deallocate(pattern2din) endif - do nm=1,isize ! loop over elements since mp_bcast cannot handle integer array - call mp_bcst(isave(nm)) ! blast out seed - enddo - do nm=1,2*ndimspec ! loop over elements since mp_bcast cannot handle integer array - call mp_bcst(pattern2d(nm)) - enddo + call mp_bcst(isave,isize) ! blast out seed + call mp_bcst(pattern2d,2*ndimspec) call random_seed(put=isave,stat=rpattern%rstate) ! subset do nn=1,len_trie_ls