From 761304138293cd63ccfd876b6cf06f0ea569dcf3 Mon Sep 17 00:00:00 2001 From: Azadeh Date: Tue, 15 Aug 2023 21:12:43 +0000 Subject: [PATCH 01/13] Thompson_MP_modifications --- src/gsi/crtm_interface.f90 | 367 +++++++++++++++++++++++++++++++- src/gsi/general_read_gfsatm.f90 | 120 +++++++++-- src/gsi/netcdfgfs_io.f90 | 28 ++- src/gsi/prt_guess.f90 | 60 +++++- 4 files changed, 544 insertions(+), 31 deletions(-) diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90 index 2305c84340..56ffd1b758 100644 --- a/src/gsi/crtm_interface.f90 +++ b/src/gsi/crtm_interface.f90 @@ -169,7 +169,8 @@ module crtm_interface real(r_kind) , save ,allocatable,dimension(:,:) :: cloudefr ! effective radius of cloud type in CRTM real(r_kind) , save ,allocatable,dimension(:,:) :: cloud_cont ! cloud content fed into CRTM real(r_kind) , save ,allocatable,dimension(:,:) :: cloud_efr ! effective radius of cloud type in CRTM - real(r_kind) , save ,allocatable,dimension(:) :: cf ! effective radius of cloud type in CRTM + real(r_kind) , save ,allocatable,dimension(:) :: cf ! effective radius of cloud type in CRTM + real(r_kind) , save ,allocatable,dimension(:) :: ni,nr ! effective radius of cloud type in CRTM real(r_kind) , save ,allocatable,dimension(:) :: hwp_guess ! column total for each hydrometeor real(r_kind) , save ,allocatable,dimension(:) :: table,table2,tablew ! GFDL saturation water vapor pressure tables real(r_kind) , save ,allocatable,dimension(:) :: des2,desw ! GFDL saturation water vapor presure @@ -435,13 +436,17 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo allocate(cloud(nsig,n_clouds_fwd)) allocate(cloudefr(nsig,n_clouds_fwd)) allocate(icloud(n_actual_clouds)) - allocate(cf(nsig)) + allocate(cf(nsig)) + allocate(ni(nsig)) + allocate(nr(nsig)) allocate(hwp_guess(n_clouds_fwd)) cloud_cont=zero cloud_efr =zero cloud =zero cloudefr =zero - cf =zero + cf =zero + ni =zero + nr =zero hwp_guess =zero call gsi_bundlegetpointer(gsi_metguess_bundle(1),cloud_names,icloud,ier) @@ -860,8 +865,11 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo endif ! regional or IGBP ! Initial GFDL saturation water vapor pressure tables + +!Option for Thompson microphysics !Azadeh + if (use_gfdl_qsat) then - if (n_actual_aerosols_wk>0 .or. n_clouds_fwd_wk>0 .and. imp_physics==11) then + if (n_actual_aerosols_wk>0 .or. n_clouds_fwd_wk>0 .and. imp_physics==11 .or. imp_physics==08) then if (mype==0) write(6,*)myname_,':initial and load GFDL saturation water vapor pressure tables' @@ -965,7 +973,9 @@ subroutine destroy_crtm if(allocated(jcloud)) deallocate(jcloud) if(allocated(cloud_cont)) deallocate(cloud_cont) if(allocated(cloud_efr)) deallocate(cloud_efr) - if(allocated(cf)) deallocate(cf) + if(allocated(cf)) deallocate(cf) + if(allocated(ni)) deallocate(ni) + if(allocated(nr)) deallocate(nr) if(allocated(hwp_guess)) deallocate(hwp_guess) if(allocated(icw)) deallocate(icw) if(allocated(lcloud4crtm_wk)) deallocate(lcloud4crtm_wk) @@ -1077,6 +1087,8 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & use obsmod, only: iadate use aeroinfo, only: nsigaerojac use chemmod, only: lread_ext_aerosol !for separate aerosol input file + use jfunc, only: jiter + use crtm_cloudcover_define, only: cloudcover_maximum_overlap, & cloudcover_random_overlap, & cloudcover_maxran_overlap, & @@ -1131,11 +1143,14 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & integer(i_kind):: i_minus, i_plus, j_minus, j_plus integer(i_kind):: itsig,itsigp,itsfc,itsfcp,itaer,itaerp integer(i_kind):: istyp00,istyp01,istyp10,istyp11 - integer(i_kind):: iqs,iozs,icfs + integer(i_kind):: iqs,iozs,icfs + integer(i_kind):: inis,inrs integer(i_kind):: error_status_clr integer(i_kind):: idx700,dprs,dprs_min integer(i_kind),dimension(8)::obs_time,anal_time integer(i_kind),dimension(msig) :: klevel + real(r_kind),dimension(nsig) :: qsat + real(r_kind),dimension(nsig) :: qcond ! ****************************** ! Constrained indexing for lai @@ -1183,7 +1198,13 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & real(r_kind),pointer,dimension(:,:,:)::aeroges_itsig =>NULL() real(r_kind),pointer,dimension(:,:,:)::aeroges_itsigp=>NULL() real(r_kind),pointer,dimension(:,:,:)::cfges_itsig =>NULL() - real(r_kind),pointer,dimension(:,:,:)::cfges_itsigp=>NULL() + real(r_kind),pointer,dimension(:,:,:)::cfges_itsigp=>NULL() + real(r_kind),pointer,dimension(:,:,:)::niges_itsig =>NULL() + real(r_kind),pointer,dimension(:,:,:)::niges_itsigp=>NULL() + real(r_kind),pointer,dimension(:,:,:)::nrges_itsig =>NULL() + real(r_kind),pointer,dimension(:,:,:)::nrges_itsigp=>NULL() + logical :: lmfdeep2 + real(r_kind) , parameter:: xrc3=200.0_r_kind logical :: sea,icmask @@ -1208,6 +1229,8 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & cloud_efr = zero cloudefr = zero cf = zero + ni = zero + nr = zero endif dx = data_s(ilat) ! grid relative latitude @@ -1354,6 +1377,15 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & icfs=istatus call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'cf',cfges_itsigp,icfs) icfs=icfs+istatus + call gsi_bundlegetpointer(gsi_metguess_bundle(itsig ),'ni',niges_itsig,istatus) + inis=istatus + call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'ni',niges_itsigp,istatus) + inis=inis+istatus + call gsi_bundlegetpointer(gsi_metguess_bundle(itsig ),'nr',nrges_itsig,istatus) + inrs=istatus + call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'nr',nrges_itsigp,istatus) + inrs=inrs+istatus + ! Space-time interpolation of temperature (h) and q fields from sigma files !$omp parallel do schedule(dynamic,1) private(k,ii,iii) @@ -1728,6 +1760,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & ges_qsat (ixp,iy ,k, itsigp)*w10+ & ges_qsat (ix ,iyp,k, itsigp)*w01+ & ges_qsat (ixp,iyp,k, itsigp)*w11)*dtsigp + rh(k) = q(k)/qsat(k) !added for cloud fraction computation c3(k)=r1000/(one-q(k)) qmix(k)=q(k)*c3(k) !convert specific humidity to mixing ratio ! Space-time interpolation of ozone(poz) @@ -1760,6 +1793,38 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & cf(k)=min(max(zero,cf(k)),one) endif ! cf + +! Space-time interpolation of cloud fraction (cf) + if (n_clouds_fwd_wk>0 .and. inis==0) then + ni(k)=((niges_itsig (ix ,iy ,k)*w00+ & + niges_itsig (ixp,iy ,k)*w10+ & + niges_itsig (ix ,iyp,k)*w01+ & + niges_itsig (ixp,iyp,k)*w11)*dtsig + & + (niges_itsigp(ix ,iy ,k)*w00+ & + niges_itsigp(ixp,iy ,k)*w10+ & + niges_itsigp(ix ,iyp,k)*w01+ & + niges_itsigp(ixp,iyp,k)*w11)*dtsigp) + +! Ensure ozone is greater than ozsmall + + ni(k)=max(zero,ni(k)) + endif ! ni + +! Space-time interpolation of cloud fraction (cf) + if (n_clouds_fwd_wk>0 .and. inrs==0) then + nr(k)=((nrges_itsig (ix ,iy ,k)*w00+ & + nrges_itsig (ixp,iy ,k)*w10+ & + nrges_itsig (ix ,iyp,k)*w01+ & + nrges_itsig (ixp,iyp,k)*w11)*dtsig + & + (nrges_itsigp(ix ,iy ,k)*w00+ & + nrges_itsigp(ixp,iy ,k)*w10+ & + nrges_itsigp(ix ,iyp,k)*w01+ & + nrges_itsigp(ixp,iyp,k)*w11)*dtsigp) + +! Ensure ozone is greater than ozsmall + + nr(k)=max(zero,nr(k)) + endif ! nr ! Quantities required for MW cloudy radiance calculations if (n_clouds_fwd_wk>0) then @@ -1820,6 +1885,15 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & end do endif + + ! Calculate Thompson effective radius for each hydrometeor + if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==08 .and. lprecip_wk) then + do ii = 1, n_clouds_fwd_wk + iii=jcloud(ii) + call calc_thompson_reff(rho_air,h,cloud(:,ii),cloud_names(iii),cloudefr(:,ii)) + end do + endif + ! Calculate GFDL cloud fraction (if no cf in metguess table) based on PDF scheme ! if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==11 .and. lcalc_gfdl_cfrac ) then if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==11 .and. lprecip_wk ) then @@ -1829,6 +1903,18 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & icfs = 0 ! load cloud fraction into CRTM endif + + ! Calculate Thompson cloud fraction !Azadeh + if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==08 .and. lprecip_wk ) then + cf_calc = zero +! sum of the cloud condensate amount for all 5 hydrometeos (ql + qi + qs + qg + qh + qr) !Azadeh + qcond(:) = cloud(:,1) + cloud(:,2) + cloud(:,3) + cloud(:,4) + cloud(:,5) + call calc_thompson_cloudfrac (nsig , lmfdeep2, xrc3, prsl, qcond(:), rh, qsat,cf_calc) + cf = cf_calc + icfs = 0 ! load cloud fraction into CRTM + endif + + ! Interpolate level pressure to observation point for top interface prsi(nsig+1)=(ges_prsi(ix ,iy ,nsig+1,itsig )*w00+ & ges_prsi(ixp,iy ,nsig+1,itsig )*w10+ & @@ -2044,11 +2130,21 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & do ii=1,n_clouds_fwd_wk !cloud_cont(k,ii)=cloud(kk2,ii)*kgkg_kgm2 cloud_cont(k,ii)=cloud(kk2,ii)*c6(k) - if (imp_physics==11 .and. lprecip_wk .and. cloud_cont(k,ii) > 1.0e-6_r_kind) then + !Option for Thompson added + if (imp_physics==11 .or. imp_physics==08 .and. lprecip_wk .and. cloud_cont(k,ii) > 1.0e-6_r_kind) then cloud_efr(k,ii)=cloudefr(kk2,ii) else cloud_efr(k,ii)=zero endif + if (trim(cloud_names_fwd(ii))=='ni' .and. atmosphere(1)%temperature(k) qmin) then + mu_w = MAX(2, MIN((NINT(1000.E6/Nt_c) + 2), 15)) + lam_w=exp(1.0_r_kind / 3.0_r_kind * log ((am_w*Nt_c *gamma(mu_w + 3.0_r_kind + 1))/(qx*gamma(mu_w+1)))) + + reff(k) = 0.5_r_kind * ((3.0_r_kind+mu_w)/lam_w)*1.0e6_r_kind + reff(k) = max(reff_min, min(reff_max, reff(k))) + else + reff(k) = zero + endif + enddo + + ! Cloud Ice + else if (trim(cloud_name)=='qi') then + am_i = rho_i*pi/6.0 + reff_min = reff_i_min + reff_max = reff_i_max + do k = 1, nsig + qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) + if (qx > qmin) then + lam_i=exp(1.0_r_kind / 3.0_r_kind * log((am_i*ni(k) *gamma(mu_i + 3.0_r_kind + 1))/(qx*gamma(mu_i+1)))) + reff(k) = 0.5_r_kind * (3.0_r_kind /lam_i)*1.0e6_r_kind + reff(k) = max(reff_min, min(reff_max, reff(k))) + else + reff(k) = zero + endif + enddo + !Rain + else if (trim(cloud_name)=='qr') then + am_r = rho_r*pi/6.0 + reff_min = reff_r_min + reff_max = reff_r_max + do k = 1, nsig + qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) + if (qx > qmin) then + lam_r=exp(1.0_r_kind / 3.0_r_kind * log ((am_r*nr(k) *gamma(mu_r + 3.0_r_kind + 1))/(qx*gamma(mu_r+1)))) + reff(k) = 0.5_r_kind *(3.0_r_kind/lam_r)*1.0e6_r_kind + reff(k) = max(reff_min, min(reff_max, reff(k))) + else + reff(k) = zero + endif + enddo + +! Snow (Field et al. 2005) + + else if (trim(cloud_name)=='qs') then + reff_min = reff_s_min + reff_max = reff_s_max + do k = 1, nsig +!..Calculate bm_s+1 (th) moment, smoc. Useful for diameter calcs. + tc0 = MIN(-0.1, tsen(k)-273.15) + qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) + if (qx > qmin) then + smob =qx/am_s + loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) & + + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 & + + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) & + + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 & + + sa(10)*cse(1)*cse(1)*cse(1) + a_ = 10.0**loga_ + b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) & + + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) & + + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) & + + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1) + smoc = a_ * smob**b_ + reff(k) = MAX(2.51E-6, MIN(0.5*(smoc/smob), 1999.E-6)) + reff(k) = max(reff_min, min(reff_max, reff(k)*1.0e6_r_kind)) + else + reff(k) = zero + endif + enddo + + ! Graupel + else if (trim(cloud_name)=='qg') then + reff_min = reff_g_min + reff_max = reff_g_max + am_g = rho_g*pi/6.0 + do k = 1, nsig + qx = qxmr(k)*rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) + if (qx > qmin) then + lam_exp = no_exp* exp(1/(3.0_r_kind+1) * log ((am_g*gamma(3.0_r_kind +1))/qx)) + lam_g=lam_exp*exp(1/3.0_r_kind*log(gamma(3.0_r_kind+mu_g+1)/((3.0_r_kind+mu_g+1)*(mu_g+1)))) + reff(k) = 3.0_r_kind/ lam_g + reff(k) = max(reff_min, min(reff_max, reff(k))) + else + reff(k) = zero + endif + enddo + + ! Mysterious + else + call die(myname_,"cannot recognize cloud name <"//trim(myname_)//">") + endif + + end subroutine calc_thompson_reff + subroutine calc_gfdl_cloudfrac(den,pt1,qv,cloud,hs,area,qsat,cfrac) !$$$ subprogram documentation block ! . . . . @@ -2796,6 +3071,82 @@ subroutine calc_gfdl_cloudfrac(den,pt1,qv,cloud,hs,area,qsat,cfrac) enddo !k-loop end subroutine calc_gfdl_cloudfrac + subroutine calc_thompson_cloudfrac(nsig, lmfdeep2, xrc3, prsl, clwf, rhly, qstl,cldtot) +!$$$ subprogram documentation block +! . . . . +! subprogram: calc_thompson_cloudfrac calculate cloud fractions from cloud using Thompson scheme +! +! prgmmr: Azadeh.Gholoubi +! +! abstract: +! +! program history log: +! +! 2023-01-24 +! +! language: f90 +! +! ==================== definition of variables ==================== ! +! input variables: ! +! prsl (msig) : model layer mean pressure in cb ! +! qstl (nsig) : layer saturate humidity in gm/gm ! +! rhly (nsig) : layer relative humidity (=qlyr/qstl) ! +! nsig : vertical layer/level dimensions ! +! lmfdeep2 : logical - true for mass flux deep convection ! +! clwf : sum of the all 5 Hydrometeors ! +! lmfdeep2 : scale-aware mass-flux deep conv scheme flag ! ! +! output variables: ! +! cldtot (nsig) - layer total cloud fraction ! +! ==================== end of description ===================== ! +!$$$ +!-------- + use constants + implicit none + +! --- inputs: + integer(i_kind), intent(in) :: nsig + real(r_kind) , dimension(nsig), intent(in) :: prsl, clwf, rhly, qstl + logical, intent(in) :: lmfdeep2 + real(r_kind) , intent(in):: xrc3 + +! --- outputs + real(r_kind) , dimension(nsig), intent(inout) :: cldtot + +! --- local variables: + real(r_kind) :: clwmin, clwm, clwt, onemrh, value, tem1, tem2 + integer(i_kind) :: k + +!> - Compute layer cloud fraction. + clwmin = 0.0 + do k = 1, nsig-1 +! converting prsl from cb to bar + clwt = 1.0e-10 * (prsl(k)*0.01) + if (clwf(k) > clwt) then + if(rhly(k) > 0.99) then + cldtot(k) = 1. + else + onemrh= max( 1.e-10, 1.0-rhly(k) ) + clwm = clwmin / max( 0.01, prsl(k)*0.01 ) + + tem1 = min(max((onemrh*qstl(k))**0.49,0.0001),1.0) + if (lmfdeep2) then + tem1 = xrc3 / tem1 + else + tem1 = 100.0 / tem1 + endif + + value = max( min( tem1*(clwf(k)-clwm), 50.0 ), 0.0 ) + tem2 = sqrt( sqrt(rhly(k)) ) + cldtot(k) = max( tem2*(1.0-exp(-value)), 0.0 ) + endif + else + cldtot(k) = 0.0 + endif + enddo + + end subroutine calc_thompson_cloudfrac +!............................................ + subroutine qs_table(n) !$$$ subprogram documentation block ! . . . . diff --git a/src/gsi/general_read_gfsatm.f90 b/src/gsi/general_read_gfsatm.f90 index ffdcd90c79..aee1531cc7 100755 --- a/src/gsi/general_read_gfsatm.f90 +++ b/src/gsi/general_read_gfsatm.f90 @@ -189,7 +189,7 @@ subroutine general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, end subroutine general_reload subroutine general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vdflag,g_cf) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vdflag,g_ni,g_nr,g_cf) ! !USES: @@ -211,7 +211,8 @@ subroutine general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & real(r_kind),dimension(grd%lat2,grd%lon2), intent( out) :: g_ps real(r_kind),dimension(grd%lat2,grd%lon2), intent(inout) :: g_z real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig),intent( out) :: g_u,g_v,& - g_vor,g_div,g_q,g_oz,g_tv,g_ql,g_qi,g_qr,g_qs,g_qg + g_vor,g_div,g_q,g_oz,g_tv,g_ql,g_qi,g_qr,g_qs,g_qg,g_ni,g_nr + real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig),intent( out),optional :: g_cf @@ -392,7 +393,25 @@ subroutine general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & g_qg(i,j,klev)=sub(ij,k) enddo enddo - elseif ( iflag(k) == 15 .and. present(g_cf) ) then + elseif ( iflag(k) == 15 ) then + klev=ilev(k) + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_ni(i,j,klev)=sub(ij,k) + enddo + enddo + elseif ( iflag(k) == 16 ) then + klev=ilev(k) + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_nr(i,j,klev)=sub(ij,k) + enddo + enddo + elseif ( iflag(k) == 17 .and. present(g_cf) ) then klev=ilev(k) ij=0 do j=1,grd%lon2 @@ -476,6 +495,7 @@ subroutine general_reload_sfc(grd,g_t2m, g_q2m,g_ps,icount,iflag,work) enddo ! do k=1,icount icount=0 + iflag=0 return @@ -2827,7 +2847,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z real(r_kind),pointer,dimension(:,:) :: g_t2m, g_q2m real(r_kind),pointer,dimension(:,:,:) :: g_vor,g_div,& g_q,g_oz,g_tv - real(r_kind),pointer,dimension(:,:,:) :: g_ql,g_qi,g_qr,g_qs,g_qg + real(r_kind),pointer,dimension(:,:,:) :: g_ql,g_qi,g_qr,g_qs,g_qg,g_ni,g_nr real(r_kind),allocatable,dimension(:,:) :: g_z real(r_kind),allocatable,dimension(:,:,:) :: g_u,g_v @@ -3037,6 +3057,8 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z call gsi_bundlegetpointer(gfs_bundle,'qs',g_qs ,ier);istatus1=istatus1+ier call gsi_bundlegetpointer(gfs_bundle,'qg',g_qg ,ier);istatus1=istatus1+ier ! call gsi_bundlegetpointer(gfs_bundle,'cf',g_cf ,ier);istatus1=istatus1+ier + call gsi_bundlegetpointer(gfs_bundle,'ni',g_ni ,ier);istatus1=istatus1+ier + call gsi_bundlegetpointer(gfs_bundle,'nr',g_nr ,ier);istatus1=istatus1+ier if ( istatus1 /= 0 ) then if ( mype == 0 ) then write(6,*) 'general_read_gfsatm_allhydro_nc: ERROR' @@ -3099,7 +3121,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif endif @@ -3134,7 +3156,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif ! Thermodynamic variable: s-->g transform, communicate to all tasks @@ -3171,7 +3193,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif end do @@ -3228,7 +3250,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif end do do k=1,nlevs @@ -3284,7 +3306,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif end do @@ -3320,7 +3342,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif icount=icount+1 @@ -3351,7 +3373,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif end do endif ! if ( uvflag ) @@ -3382,7 +3404,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif end do @@ -3413,7 +3435,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif end do @@ -3444,7 +3466,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif enddo ! do k=1,nlevs @@ -3474,7 +3496,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif enddo ! do k=1,nlevs @@ -3504,7 +3526,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif enddo ! do k=1,nlevs @@ -3534,7 +3556,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif enddo ! do k=1,nlevs @@ -3564,13 +3586,73 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm .or. k==nlevs) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) + endif + enddo ! do k=1,nlevs + + do k=1,nlevs + icount=icount+1 + iflag(icount)=15 + ilev(icount)=k + kr = levs+1-k ! netcdf is top to bottom, need to flip + + if (mype==mype_use(icount)) then + call read_vardata(filges, 'nccice', rwork3d0, nslice=kr, slicedim=3) + ! cloud ice water number concentration. + if ( diff_res ) then + grid_b=rwork3d0(:,:,1) + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork3d0(:,:,1) + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm .or. k==nlevs ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) + endif + enddo ! do k=1,nlevs + + do k=1,nlevs + icount=icount+1 + iflag(icount)=16 + ilev(icount)=k + kr = levs+1-k ! netcdf is top to bottom, need to flip + + if (mype==mype_use(icount)) then + call read_vardata(filges, 'nconrd', rwork3d0, nslice=kr, slicedim=3) + ! rain number concentration. + if ( diff_res ) then + grid_b=rwork3d0(:,:,1) + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork3d0(:,:,1) + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm .or. k==nlevs) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif enddo ! do k=1,nlevs ! do k=1,nlevs ! icount=icount+1 -! iflag(icount)=15 +! iflag(icount)=17 ! ilev(icount)=k ! kr = levs+1-k ! netcdf is top to bottom, need to flip ! diff --git a/src/gsi/netcdfgfs_io.f90 b/src/gsi/netcdfgfs_io.f90 index 41e8f33e03..be0c1c1cce 100644 --- a/src/gsi/netcdfgfs_io.f90 +++ b/src/gsi/netcdfgfs_io.f90 @@ -159,6 +159,9 @@ subroutine read_ real(r_kind),pointer,dimension(:,:,:):: ges_qs_it => NULL() real(r_kind),pointer,dimension(:,:,:):: ges_qg_it => NULL() real(r_kind),pointer,dimension(:,:,:):: ges_cf_it => NULL() + real(r_kind),pointer,dimension(:,:,:):: ges_ni_it => NULL() + real(r_kind),pointer,dimension(:,:,:):: ges_nr_it => NULL() + type(sub2grid_info) :: grd_t logical regional @@ -169,7 +172,7 @@ subroutine read_ integer(i_kind),parameter :: n2d=2 ! integer(i_kind),parameter :: n3d=8 integer(i_kind),parameter :: n2d_2m=4 - integer(i_kind),parameter :: n3d=14 + integer(i_kind),parameter :: n3d=16 character(len=4), parameter :: vars2d(n2d) = (/ 'z ', 'ps ' /) character(len=4), parameter :: vars2d_with2m(n2d_2m) = (/ 'z ', 'ps ','t2m ','q2m ' /) ! character(len=4), parameter :: vars3d(n3d) = (/ 'u ', 'v ', & @@ -182,13 +185,17 @@ subroutine read_ 'cw ', 'oz ', & 'ql ', 'qi ', & 'qr ', 'qs ', & - 'qg ', 'cf ' /) + 'qg ', 'ni ', & + 'nr ', 'cf ' /) + real(r_kind),pointer,dimension(:,:):: ptr2d =>NULL() real(r_kind),pointer,dimension(:,:,:):: ptr3d =>NULL() regional=.false. inner_vars=1 - num_fields=min(14*grd_a%nsig+2,npe) + + num_fields=min(n3d*grd_a%nsig+2,npe) +! num_fields=min(14*grd_a%nsig+2,npe) ! Create temporary communication information fore read routines call general_sub2grid_create_info(grd_t,inner_vars,grd_a%nlat,grd_a%nlon, & grd_a%nsig,num_fields,regional) @@ -200,6 +207,7 @@ subroutine read_ else call gsi_bundlecreate(atm_bundle,atm_grid,'aux-atm-read',istatus,names2d=vars2d,names3d=vars3d) endif + if(istatus/=0) then write(6,*) myname_,': trouble creating atm_bundle' call stop2(999) @@ -248,6 +256,8 @@ subroutine read_ if (associated(ges_qr_it)) ges_qr_it(i,j,k) = max(qcmin,ges_qr_it(i,j,k)) if (associated(ges_qs_it)) ges_qs_it(i,j,k) = max(qcmin,ges_qs_it(i,j,k)) if (associated(ges_qg_it)) ges_qg_it(i,j,k) = max(qcmin,ges_qg_it(i,j,k)) + if (associated(ges_ni_it)) ges_ni_it(i,j,k) = max(qcmin,ges_ni_it(i,j,k)) + if (associated(ges_nr_it)) ges_nr_it(i,j,k) = max(qcmin,ges_nr_it(i,j,k)) if (associated(ges_cf_it)) ges_cf_it(i,j,k) = min(max(zero,ges_cf_it(i,j,k)),one) enddo enddo @@ -256,6 +266,8 @@ subroutine read_ l_cld_derived = associated(ges_cwmr_it).and.& associated(ges_q_it) .and.& associated(ges_ql_it) .and.& + associated(ges_ni_it) .and.& + associated(ges_nr_it) .and.& associated(ges_qi_it) .and.& associated(ges_tv_it) ! call set_cloud_lower_bound(ges_cwmr_it) @@ -343,6 +355,16 @@ subroutine set_guess_ call gsi_bundlegetpointer (gsi_metguess_bundle(it),'ql',ges_ql_it,istatus) if(istatus==0) ges_ql_it = ptr3d endif + call gsi_bundlegetpointer (atm_bundle,'ni',ptr3d,istatus) + if (istatus==0) then + call gsi_bundlegetpointer (gsi_metguess_bundle(it),'ni',ges_ni_it,istatus) + if(istatus==0) ges_ni_it = ptr3d + endif + call gsi_bundlegetpointer (atm_bundle,'nr',ptr3d,istatus) + if (istatus==0) then + call gsi_bundlegetpointer (gsi_metguess_bundle(it),'nr',ges_nr_it,istatus) + if(istatus==0) ges_nr_it = ptr3d + endif call gsi_bundlegetpointer (atm_bundle,'qi',ptr3d,istatus) if (istatus==0) then call gsi_bundlegetpointer (gsi_metguess_bundle(it),'qi',ges_qi_it,istatus) diff --git a/src/gsi/prt_guess.f90 b/src/gsi/prt_guess.f90 index 97c1e0a7c1..11befcd2da 100644 --- a/src/gsi/prt_guess.f90 +++ b/src/gsi/prt_guess.f90 @@ -287,7 +287,8 @@ subroutine prt_guess2(sgrep) integer(i_kind), parameter :: nvars3=5 integer(i_kind) :: nvars,nvars2,nvarsc,nc integer(i_kind) ii,istatus,ier,ivar - integer(i_kind) iql,iqi,iqr,iqs,iqg,icf + integer(i_kind) iql,iqi,iqr,iqs,iqg,icf + integer(i_kind) ini,inr integer(i_kind) ntsig integer(i_kind) ntsfc integer(i_kind) n_actual_clouds @@ -307,6 +308,8 @@ subroutine prt_guess2(sgrep) real(r_kind),pointer,dimension(:,:,:)::ges_qr_it => NULL() real(r_kind),pointer,dimension(:,:,:)::ges_qs_it => NULL() real(r_kind),pointer,dimension(:,:,:)::ges_qg_it => NULL() + real(r_kind),pointer,dimension(:,:,:)::ges_ni_it => NULL() + real(r_kind),pointer,dimension(:,:,:)::ges_nr_it => NULL() real(r_kind),pointer,dimension(:,:,:)::ges_cf_it => NULL() ! character(len=4) :: cvar(nvars+3) character(len=4),allocatable,dimension(:) :: cvar @@ -345,26 +348,43 @@ subroutine prt_guess2(sgrep) if (ivar > 0) then call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'ql',ges_ql_it,istatus) ier=ier+istatus + if (ier==0) nvarsc=nvarsc+1 endif call gsi_metguess_get ( 'var::qi', ivar, ier ); iqi=ivar if (ivar > 0) then call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'qi',ges_qi_it,istatus) ier=ier+istatus + if (ier==0) nvarsc=nvarsc+1 endif call gsi_metguess_get ( 'var::qr', ivar, ier ); iqr=ivar if (ivar > 0) then call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'qr',ges_qr_it,istatus) ier=ier+istatus + if (ier==0) nvarsc=nvarsc+1 endif call gsi_metguess_get ( 'var::qs', ivar, ier ); iqs=ivar if (ivar > 0) then call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'qs',ges_qs_it,istatus) ier=ier+istatus + if (ier==0) nvarsc=nvarsc+1 endif call gsi_metguess_get ( 'var::qg', ivar, ier ); iqg=ivar if (ivar > 0) then call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'qg',ges_qg_it,istatus) ier=ier+istatus + if (ier==0) nvarsc=nvarsc+1 + endif + call gsi_metguess_get ( 'var::ni', ivar, ier ); ini=ivar + if (ivar > 0) then + call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'ni',ges_ni_it,istatus) + ier=ier+istatus + if (ier==0) nvarsc=nvarsc+1 + endif + call gsi_metguess_get ( 'var::nr', ivar, ier ); inr=ivar + if (ivar > 0) then + call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'nr',ges_nr_it,istatus) + ier=ier+istatus + if (ier==0) nvarsc=nvarsc+1 endif end if nvarsc=n_actual_clouds @@ -415,6 +435,14 @@ subroutine prt_guess2(sgrep) nc=nc+1 cvar(nvars1+nc)='QG ' endif + if(ini>0) then + nc=nc+1 + cvar(nvars1+nc)='NI ' + endif + if(inr>0) then + nc=nc+1 + cvar(nvars1+nc)='NR ' + endif if(icf>0) then nc=nc+1 cvar(nvars1+nc)='CF ' @@ -455,6 +483,17 @@ subroutine prt_guess2(sgrep) nc=nc+1 zloc(nvars1+nc) = sum (ges_qg_it (2:lat1+1,2:lon1+1,1:nsig)) endif + + if(ini>0) then + nc=nc+1 + zloc(nvars1+nc) = sum (ges_ni_it (2:lat1+1,2:lon1+1,1:nsig)) + endif + + if(inr>0) then + nc=nc+1 + zloc(nvars1+nc) = sum (ges_nr_it (2:lat1+1,2:lon1+1,1:nsig)) + endif + if(icf>0) then nc=nc+1 zloc(nvars1+nc) = sum (ges_cf_it (2:lat1+1,2:lon1+1,1:nsig)) @@ -491,6 +530,16 @@ subroutine prt_guess2(sgrep) nc=nc+1 zloc(nvars+nvars1+nc) = minval(ges_qg_it (2:lat1+1,2:lon1+1,1:nsig)) endif + + if(ini>0) then + nc=nc+1 + zloc(nvars+nvars1+nc) = minval(ges_ni_it (2:lat1+1,2:lon1+1,1:nsig)) + endif + if(inr>0) then + nc=nc+1 + zloc(nvars+nvars1+nc) = minval(ges_nr_it (2:lat1+1,2:lon1+1,1:nsig)) + endif + if(icf>0) then nc=nc+1 zloc(nvars+nvars1+nc) = minval(ges_cf_it (2:lat1+1,2:lon1+1,1:nsig)) @@ -527,6 +576,15 @@ subroutine prt_guess2(sgrep) nc=nc+1 zloc(2*nvars+nvars1+nc) = maxval(ges_qg_it (2:lat1+1,2:lon1+1,1:nsig)) endif + if(ini>0) then + nc=nc+1 + zloc(2*nvars+nvars1+nc) = maxval(ges_ni_it (2:lat1+1,2:lon1+1,1:nsig)) + endif + if(inr>0) then + nc=nc+1 + zloc(2*nvars+nvars1+nc) = maxval(ges_nr_it (2:lat1+1,2:lon1+1,1:nsig)) + endif + if(icf>0) then nc=nc+1 zloc(2*nvars+nvars1+nc) = maxval(ges_cf_it (2:lat1+1,2:lon1+1,1:nsig)) From 7c6bc78fb8ba557b961ab99a0c721661be1a9b7f Mon Sep 17 00:00:00 2001 From: Azadeh Date: Thu, 4 Apr 2024 18:43:47 +0000 Subject: [PATCH 02/13] Modified code based on FV3jedi code --- src/gsi/crtm_interface.f90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90 index 3a08ba7260..ea39b01e68 100644 --- a/src/gsi/crtm_interface.f90 +++ b/src/gsi/crtm_interface.f90 @@ -2748,8 +2748,8 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) real(r_kind), parameter :: reff_r_max = 1000.0_r_kind ! [micron] ! previous value was 10000.0_r_kind real(r_kind), parameter:: mu_r = 0.0 ! Parameters for snow - real(r_kind), parameter :: reff_s_min = 50.0_r_kind ! [micron] ! previous value was 0.0_r_kind - real(r_kind), parameter :: reff_s_max = 1000.0_r_kind ! [micron] ! previous value was 10000.0_r_kind + real(r_kind), parameter :: reff_s_min = 5.0_r_kind ! [micron] ! previous value was 0.0_r_kind + real(r_kind), parameter :: reff_s_max = 5000.0_r_kind ! [micron] ! previous value was 10000.0_r_kind !For snow moments conversions (from Field et al. 2005) real, dimension(10), parameter:: & @@ -2777,8 +2777,8 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) reff_max = reff_w_max do k = 1, nsig qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) - if (qx > qmin) then - mu_w = MAX(2, MIN((NINT(1000.E6/Nt_c) + 2), 15)) + if (qx < qmin) then + mu_w = MAX(2, MIN((NINT(1000.E6/Nt_c) + 2), 15)) lam_w=exp(1.0_r_kind / 3.0_r_kind * log ((am_w*Nt_c *gamma(mu_w + 3.0_r_kind + 1))/(qx*gamma(mu_w+1)))) reff(k) = 0.5_r_kind * ((3.0_r_kind+mu_w)/lam_w)*1.0e6_r_kind @@ -2795,7 +2795,7 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) reff_max = reff_i_max do k = 1, nsig qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) - if (qx > qmin) then + if (qx < qmin) then lam_i=exp(1.0_r_kind / 3.0_r_kind * log((am_i*ni(k) *gamma(mu_i + 3.0_r_kind + 1))/(qx*gamma(mu_i+1)))) reff(k) = 0.5_r_kind * (3.0_r_kind /lam_i)*1.0e6_r_kind reff(k) = max(reff_min, min(reff_max, reff(k))) @@ -2810,7 +2810,7 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) reff_max = reff_r_max do k = 1, nsig qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) - if (qx > qmin) then + if (qx < qmin) then lam_r=exp(1.0_r_kind / 3.0_r_kind * log ((am_r*nr(k) *gamma(mu_r + 3.0_r_kind + 1))/(qx*gamma(mu_r+1)))) reff(k) = 0.5_r_kind *(3.0_r_kind/lam_r)*1.0e6_r_kind reff(k) = max(reff_min, min(reff_max, reff(k))) @@ -2828,7 +2828,7 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) !..Calculate bm_s+1 (th) moment, smoc. Useful for diameter calcs. tc0 = MIN(-0.1, tsen(k)-273.15) qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) - if (qx > qmin) then + if (qx < qmin) then smob =qx/am_s loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) & + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 & @@ -2855,7 +2855,7 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) am_g = rho_g*pi/6.0 do k = 1, nsig qx = qxmr(k)*rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) - if (qx > qmin) then + if (qx < qmin) then lam_exp = no_exp* exp(1/(3.0_r_kind+1) * log ((am_g*gamma(3.0_r_kind +1))/qx)) lam_g=lam_exp*exp(1/3.0_r_kind*log(gamma(3.0_r_kind+mu_g+1)/((3.0_r_kind+mu_g+1)*(mu_g+1)))) reff(k) = 3.0_r_kind/ lam_g From 3a197ca5521ec88e71ad37b75b4c5868300d4d69 Mon Sep 17 00:00:00 2001 From: Azadeh Date: Fri, 5 Apr 2024 01:15:19 +0000 Subject: [PATCH 03/13] Corrected g ease enter the commit message for your changes. Lines starting --- src/gsi/crtm_interface.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90 index ea39b01e68..dc4a0a944a 100644 --- a/src/gsi/crtm_interface.f90 +++ b/src/gsi/crtm_interface.f90 @@ -2857,8 +2857,8 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) qx = qxmr(k)*rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) if (qx < qmin) then lam_exp = no_exp* exp(1/(3.0_r_kind+1) * log ((am_g*gamma(3.0_r_kind +1))/qx)) - lam_g=lam_exp*exp(1/3.0_r_kind*log(gamma(3.0_r_kind+mu_g+1)/((3.0_r_kind+mu_g+1)*(mu_g+1)))) - reff(k) = 3.0_r_kind/ lam_g + lam_g=lam_exp*exp(1/3.0_r_kind*log(gamma(3.0_r_kind+mu_g+1)/((3.0_r_kind+mu_g+1)*(mu_g+1)))) + reff(k) = 0.5_r_kind *(3.0_r_kind/ lam_g)*1.0e6_r_kind reff(k) = max(reff_min, min(reff_max, reff(k))) else reff(k) = zero From 7623bac321e419b296ac0a3b8950225e12a73485 Mon Sep 17 00:00:00 2001 From: Azadeh Date: Wed, 8 May 2024 19:53:17 +0000 Subject: [PATCH 04/13] Addressed comments, Followed GSI code Standards --- src/gsi/crtm_interface.f90 | 109 +++++++++++++++----------------- src/gsi/general_read_gfsatm.f90 | 1 - src/gsi/netcdfgfs_io.f90 | 4 -- 3 files changed, 50 insertions(+), 64 deletions(-) diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90 index dc4a0a944a..d095c212c9 100644 --- a/src/gsi/crtm_interface.f90 +++ b/src/gsi/crtm_interface.f90 @@ -169,8 +169,8 @@ module crtm_interface real(r_kind) , save ,allocatable,dimension(:,:) :: cloudefr ! effective radius of cloud type in CRTM real(r_kind) , save ,allocatable,dimension(:,:) :: cloud_cont ! cloud content fed into CRTM real(r_kind) , save ,allocatable,dimension(:,:) :: cloud_efr ! effective radius of cloud type in CRTM - real(r_kind) , save ,allocatable,dimension(:) :: cf ! effective radius of cloud type in CRTM - real(r_kind) , save ,allocatable,dimension(:) :: ni,nr ! effective radius of cloud type in CRTM + real(r_kind) , save ,allocatable,dimension(:) :: cf ! cloud fraction + real(r_kind) , save ,allocatable,dimension(:) :: ni,nr ! ice and rain number concentration real(r_kind) , save ,allocatable,dimension(:) :: hwp_guess ! column total for each hydrometeor real(r_kind) , save ,allocatable,dimension(:) :: table,table2,tablew ! GFDL saturation water vapor pressure tables real(r_kind) , save ,allocatable,dimension(:) :: des2,desw ! GFDL saturation water vapor presure @@ -866,7 +866,7 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo ! Initial GFDL saturation water vapor pressure tables -!Option for Thompson microphysics !Azadeh +!Option for Thompson microphysics if (use_gfdl_qsat) then if (n_actual_aerosols_wk>0 .or. n_clouds_fwd_wk>0 .and. imp_physics==11 .or. imp_physics==08) then @@ -1086,8 +1086,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & use crtm_module, only: limit_exp,o3_id,toa_pressure use obsmod, only: iadate use aeroinfo, only: nsigaerojac - use chemmod, only: lread_ext_aerosol !for separate aerosol input file - use jfunc, only: jiter + use chemmod, only: lread_ext_aerosol !for separate aerosol input file use crtm_cloudcover_define, only: cloudcover_maximum_overlap, & cloudcover_random_overlap, & @@ -1762,7 +1761,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & ges_qsat (ixp,iy ,k, itsigp)*w10+ & ges_qsat (ix ,iyp,k, itsigp)*w01+ & ges_qsat (ixp,iyp,k, itsigp)*w11)*dtsigp - rh(k) = q(k)/qsat(k) !added for cloud fraction computation + rh(k) = q(k)/qsat(k) !added for cloud fraction computation c3(k)=r1000/(one-q(k)) qmix(k)=q(k)*c3(k) !convert specific humidity to mixing ratio ! Space-time interpolation of ozone(poz) @@ -1796,7 +1795,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & cf(k)=min(max(zero,cf(k)),one) endif ! cf -! Space-time interpolation of cloud fraction (cf) +! Space-time interpolation of ice number concentration (ni) if (n_clouds_fwd_wk>0 .and. inis==0) then ni(k)=((niges_itsig (ix ,iy ,k)*w00+ & niges_itsig (ixp,iy ,k)*w10+ & @@ -1806,13 +1805,10 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & niges_itsigp(ixp,iy ,k)*w10+ & niges_itsigp(ix ,iyp,k)*w01+ & niges_itsigp(ixp,iyp,k)*w11)*dtsigp) - -! Ensure ozone is greater than ozsmall - ni(k)=max(zero,ni(k)) endif ! ni -! Space-time interpolation of cloud fraction (cf) +! Space-time interpolation of rain number concentration (nr) if (n_clouds_fwd_wk>0 .and. inrs==0) then nr(k)=((nrges_itsig (ix ,iy ,k)*w00+ & nrges_itsig (ixp,iy ,k)*w10+ & @@ -1822,9 +1818,6 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & nrges_itsigp(ixp,iy ,k)*w10+ & nrges_itsigp(ix ,iyp,k)*w01+ & nrges_itsigp(ixp,iyp,k)*w11)*dtsigp) - -! Ensure ozone is greater than ozsmall - nr(k)=max(zero,nr(k)) endif ! nr ! Quantities required for MW cloudy radiance calculations @@ -1892,12 +1885,12 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==08 .and. lprecip_wk) then do ii = 1, n_clouds_fwd_wk iii=jcloud(ii) - call calc_thompson_reff(rho_air,h,cloud(:,ii),cloud_names(iii),cloudefr(:,ii)) + call calc_thompson_reff(rho_air,h,cloud(:,ii),cloud_names(iii),cloudefr(:,ii)) end do endif ! Calculate GFDL cloud fraction (if no cf in metguess table) based on PDF scheme -! if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==11 .and. lcalc_gfdl_cfrac ) then + ! if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==11 .and. lcalc_gfdl_cfrac ) then if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==11 .and. lprecip_wk ) then cf_calc = zero call calc_gfdl_cloudfrac(rho_air,h,q,cloud,hs,garea,qs,cf_calc) @@ -1905,18 +1898,15 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & icfs = 0 ! load cloud fraction into CRTM endif - - ! Calculate Thompson cloud fraction !Azadeh + ! Calculate Thompson cloud fraction if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==08 .and. lprecip_wk ) then cf_calc = zero -! sum of the cloud condensate amount for all 5 hydrometeos (ql + qi + qs + qg + qh + qr) !Azadeh + ! sum of the cloud condensate amount for all 5 hydrometeos (ql + qi + qs + qg + qh + qr) qcond(:) = cloud(:,1) + cloud(:,2) + cloud(:,3) + cloud(:,4) + cloud(:,5) - call calc_thompson_cloudfrac (nsig , lmfdeep2, xrc3, prsl, qcond(:), rh, qsat,cf_calc) + call calc_thompson_cloudfrac (nsig , lmfdeep2, xrc3, prsl, qcond(:), rh, qsat,cf_calc) cf = cf_calc icfs = 0 ! load cloud fraction into CRTM endif - - ! Interpolate level pressure to observation point for top interface prsi(nsig+1)=(ges_prsi(ix ,iy ,nsig+1,itsig )*w00+ & ges_prsi(ixp,iy ,nsig+1,itsig )*w10+ & @@ -2721,14 +2711,15 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) character(len=*), parameter :: myname_ = 'calc_thompson_reff' integer(i_kind) :: k integer(i_kind) :: j + integer(i_kind) :: mu_w real(r_kind) :: qx real(r_kind) :: reff_min, reff_max - real(r_kind):: lam_i,lam_w,lam_s,lam_r, lam_g,lam_exp, am_r,am_w,am_i,am_g, N_0r ,N_0i,N_0w,mu_w - + real(r_kind):: lam_i,lam_w,lam_s,lam_r, lam_g,lam_exp, am_r,am_w,am_i,am_g, N_0r ,N_0i,N_0w + ! Parameters real(r_kind), parameter :: qmin = 1.0e-12_r_kind ! [kg/kg ] !.. droplet number concentration. - real, parameter :: Nt_c = 100.E6 ![m-3] + real, parameter :: Nt_c = 100.E6_r_kind ![m-3] ! Parameters for water cloud real(r_kind), parameter :: ccn = 1.0e8_r_kind @@ -2740,46 +2731,46 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) real(r_kind), parameter :: rho_i = 890.0_r_kind ! [kg/m3 ] real(r_kind), parameter :: reff_i_min = 2.5_r_kind ! previous value was 10_r_kind real(r_kind), parameter :: reff_i_max = 250.0_r_kind ! previous value was 150_r_kind - real(r_kind), parameter:: mu_i = 0.0 + real(r_kind), parameter:: mu_i = 0.0_r_kind ! Parameters for rain (Lin 1983) real(r_kind), parameter :: rho_r = 1000.0_r_kind ! [kg/m3 ] real(r_kind), parameter :: reff_r_min = 50.0_r_kind ! [micron] ! previous value was 0.0_r_kind real(r_kind), parameter :: reff_r_max = 1000.0_r_kind ! [micron] ! previous value was 10000.0_r_kind - real(r_kind), parameter:: mu_r = 0.0 + real(r_kind), parameter:: mu_r = 0.0_r_kind ! Parameters for snow real(r_kind), parameter :: reff_s_min = 5.0_r_kind ! [micron] ! previous value was 0.0_r_kind real(r_kind), parameter :: reff_s_max = 5000.0_r_kind ! [micron] ! previous value was 10000.0_r_kind !For snow moments conversions (from Field et al. 2005) - real, dimension(10), parameter:: & - sa = (/ 5.065339, -0.062659, -3.032362, 0.029469, -0.000285, & - & 0.31255, 0.000204, 0.003199, 0.0, -0.015952/) - real, dimension(10), parameter:: & - sb = (/ 0.476221, -0.015896, 0.165977, 0.007468, -0.000141, & - & 0.060366, 0.000079, 0.000594, 0.0, -0.003577/) + real(r_kind), dimension(10), parameter:: & + sa = (/ 5.065339_r_kind, -0.062659_r_kind, -3.032362_r_kind, 0.029469_r_kind, -0.000285_r_kind, & + & 0.31255_r_kind, 0.000204_r_kind, 0.003199_r_kind, 0.0_r_kind, -0.015952_r_kind/) + real(r_kind), dimension(10), parameter:: & + sb = (/ 0.476221_r_kind, -0.015896_r_kind, 0.165977_r_kind, 0.007468_r_kind, -0.000141_r_kind, & + & 0.060366_r_kind, 0.000079_r_kind, 0.000594_r_kind, 0.0_r_kind, -0.003577_r_kind/) real(r_kind), parameter :: am_s = 0.069_r_kind real(r_kind), parameter :: bm_s = 2.0_r_kind - real, dimension(1), parameter :: cse = (/ bm_s + 1. /) - real :: tc0, smo2, smob, smoc, smoz, a_, b_, loga_ + real(r_kind), dimension(1), parameter :: cse = (/ bm_s + 1.0_r_kind /) + real(r_kind) :: tc0, smo2, smob, smoc, smoz, a_, b_, loga_ ! Parameters for graupel (Lin 1983) real(r_kind), parameter :: rho_g = 500.0_r_kind ! [kg/m3 ] real(r_kind), parameter :: reff_g_min = 150.0_r_kind ! [micron] ! previous value was 0.0_r_kind real(r_kind), parameter :: reff_g_max = 5000.0_r_kind ! [micron] ! previous value was 10000.0_r_kind - real(r_kind), parameter:: mu_g = 0.0 + real(r_kind), parameter:: mu_g = 0.0_r_kind real(r_kind), parameter :: no_exp = 1.0e-4_r_kind !cloud water if (trim(cloud_name)=='ql') then - am_w = rho_w*pi/6.0 + am_w = rho_w*pi/6.0_r_kind reff_min = reff_w_min reff_max = reff_w_max do k = 1, nsig qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) if (qx < qmin) then - mu_w = MAX(2, MIN((NINT(1000.E6/Nt_c) + 2), 15)) - lam_w=exp(1.0_r_kind / 3.0_r_kind * log ((am_w*Nt_c *gamma(mu_w + 3.0_r_kind + 1))/(qx*gamma(mu_w+1)))) + mu_w = MAX(2, MIN((NINT(1000.E6_r_kind/Nt_c) + 2), 15)) + lam_w=exp(1.0_r_kind / 3.0_r_kind * log ((am_w*Nt_c *gamma(mu_w + 3.0_r_kind + 1.0_r_kind))/(qx*gamma(mu_w+1.0_r_kind)))) reff(k) = 0.5_r_kind * ((3.0_r_kind+mu_w)/lam_w)*1.0e6_r_kind reff(k) = max(reff_min, min(reff_max, reff(k))) @@ -2790,13 +2781,13 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) ! Cloud Ice else if (trim(cloud_name)=='qi') then - am_i = rho_i*pi/6.0 + am_i = rho_i*pi/6.0_r_kind reff_min = reff_i_min reff_max = reff_i_max do k = 1, nsig qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) if (qx < qmin) then - lam_i=exp(1.0_r_kind / 3.0_r_kind * log((am_i*ni(k) *gamma(mu_i + 3.0_r_kind + 1))/(qx*gamma(mu_i+1)))) + lam_i=exp(1.0_r_kind / 3.0_r_kind * log((am_i*ni(k) *gamma(mu_i + 3.0_r_kind + 1.0_r_kind))/(qx*gamma(mu_i+1.0_r_kind)))) reff(k) = 0.5_r_kind * (3.0_r_kind /lam_i)*1.0e6_r_kind reff(k) = max(reff_min, min(reff_max, reff(k))) else @@ -2805,13 +2796,13 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) enddo !Rain else if (trim(cloud_name)=='qr') then - am_r = rho_r*pi/6.0 + am_r = rho_r*pi/6.0_r_kind reff_min = reff_r_min reff_max = reff_r_max do k = 1, nsig qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) if (qx < qmin) then - lam_r=exp(1.0_r_kind / 3.0_r_kind * log ((am_r*nr(k) *gamma(mu_r + 3.0_r_kind + 1))/(qx*gamma(mu_r+1)))) + lam_r=exp(1.0_r_kind / 3.0_r_kind * log ((am_r*nr(k) *gamma(mu_r + 3.0_r_kind + 1.0_r_kind))/(qx*gamma(mu_r + 1.0_r_kind)))) reff(k) = 0.5_r_kind *(3.0_r_kind/lam_r)*1.0e6_r_kind reff(k) = max(reff_min, min(reff_max, reff(k))) else @@ -2826,7 +2817,7 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) reff_max = reff_s_max do k = 1, nsig !..Calculate bm_s+1 (th) moment, smoc. Useful for diameter calcs. - tc0 = MIN(-0.1, tsen(k)-273.15) + tc0 = MIN(-0.1_r_kind, tsen(k)-273.15_r_kind) qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) if (qx < qmin) then smob =qx/am_s @@ -2841,7 +2832,7 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1) smoc = a_ * smob**b_ - reff(k) = MAX(2.51E-6, MIN(0.5*(smoc/smob), 1999.E-6)) + reff(k) = MAX(2.51E-6_r_kind, MIN(0.5*(smoc/smob), 1999.E-6_r_kind)) reff(k) = max(reff_min, min(reff_max, reff(k)*1.0e6_r_kind)) else reff(k) = zero @@ -2852,12 +2843,12 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) else if (trim(cloud_name)=='qg') then reff_min = reff_g_min reff_max = reff_g_max - am_g = rho_g*pi/6.0 + am_g = rho_g*pi/6.0_r_kind do k = 1, nsig qx = qxmr(k)*rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) if (qx < qmin) then - lam_exp = no_exp* exp(1/(3.0_r_kind+1) * log ((am_g*gamma(3.0_r_kind +1))/qx)) - lam_g=lam_exp*exp(1/3.0_r_kind*log(gamma(3.0_r_kind+mu_g+1)/((3.0_r_kind+mu_g+1)*(mu_g+1)))) + lam_exp = no_exp* exp(1.0_r_kind/(3.0_r_kind+1.0_r_kind) * log ((am_g*gamma(3.0_r_kind +1.0_r_kind))/qx)) + lam_g=lam_exp*exp(1/3.0_r_kind*log(gamma(3.0_r_kind+mu_g+1.0_r_kind)/((3.0_r_kind+mu_g+1.0_r_kind)*(mu_g+1.0_r_kind)))) reff(k) = 0.5_r_kind *(3.0_r_kind/ lam_g)*1.0e6_r_kind reff(k) = max(reff_min, min(reff_max, reff(k))) else @@ -3131,30 +3122,30 @@ subroutine calc_thompson_cloudfrac(nsig, lmfdeep2, xrc3, prsl, clwf, rhly, qstl, integer(i_kind) :: k !> - Compute layer cloud fraction. - clwmin = 0.0 + clwmin = 0.0_r_kind do k = 1, nsig-1 ! converting prsl from cb to bar - clwt = 1.0e-10 * (prsl(k)*0.01) + clwt = 1.0e-10_r_kind * (prsl(k)*0.01_r_kind) if (clwf(k) > clwt) then - if(rhly(k) > 0.99) then - cldtot(k) = 1. + if(rhly(k) > 0.99_r_kind) then + cldtot(k) = 1.0_r_kind else - onemrh= max( 1.e-10, 1.0-rhly(k) ) - clwm = clwmin / max( 0.01, prsl(k)*0.01 ) + onemrh= max( 1.e-10_r_kind, 1.0_r_kind-rhly(k) ) + clwm = clwmin / max( 0.01_r_kind, prsl(k)*0.01_r_kind ) - tem1 = min(max((onemrh*qstl(k))**0.49,0.0001),1.0) + tem1 = min(max((onemrh*qstl(k))**0.49_r_kind,0.0001_r_kind),1.0_r_kind) if (lmfdeep2) then tem1 = xrc3 / tem1 else - tem1 = 100.0 / tem1 + tem1 = 100.0_r_kind / tem1 endif - value = max( min( tem1*(clwf(k)-clwm), 50.0 ), 0.0 ) + value = max( min( tem1*(clwf(k)-clwm), 50.0_r_kind ), 0.0_r_kind ) tem2 = sqrt( sqrt(rhly(k)) ) - cldtot(k) = max( tem2*(1.0-exp(-value)), 0.0 ) + cldtot(k) = max( tem2*(1.0_r_kind-exp(-value)), 0.0_r_kind ) endif else - cldtot(k) = 0.0 + cldtot(k) = 0.0_r_kind endif enddo diff --git a/src/gsi/general_read_gfsatm.f90 b/src/gsi/general_read_gfsatm.f90 index b3f7e62b9f..1eb5ba0dd9 100755 --- a/src/gsi/general_read_gfsatm.f90 +++ b/src/gsi/general_read_gfsatm.f90 @@ -495,7 +495,6 @@ subroutine general_reload_sfc(grd,g_t2m, g_q2m,g_ps,icount,iflag,work) enddo ! do k=1,icount icount=0 - iflag=0 return diff --git a/src/gsi/netcdfgfs_io.f90 b/src/gsi/netcdfgfs_io.f90 index 57292db820..7db02636cd 100644 --- a/src/gsi/netcdfgfs_io.f90 +++ b/src/gsi/netcdfgfs_io.f90 @@ -161,12 +161,9 @@ subroutine read_ real(r_kind),pointer,dimension(:,:,:):: ges_cf_it => NULL() real(r_kind),pointer,dimension(:,:,:):: ges_ni_it => NULL() real(r_kind),pointer,dimension(:,:,:):: ges_nr_it => NULL() - - type(sub2grid_info) :: grd_t logical regional logical:: l_cld_derived,zflag,inithead - type(gsi_bundle) :: atm_bundle type(gsi_grid) :: atm_grid integer(i_kind),parameter :: n2d=2 @@ -195,7 +192,6 @@ subroutine read_ inner_vars=1 num_fields=min(n3d*grd_a%nsig+2,npe) -! num_fields=min(14*grd_a%nsig+2,npe) ! Create temporary communication information fore read routines call general_sub2grid_create_info(grd_t,inner_vars,grd_a%nlat,grd_a%nlon, & grd_a%nsig,num_fields,regional) From c24e56878aa785becae094712260ef46dc10ccc7 Mon Sep 17 00:00:00 2001 From: Azadeh Date: Thu, 9 May 2024 15:16:51 +0000 Subject: [PATCH 05/13] Updated prt_guess.f90 to print ni and nr in stdout --- src/gsi/prt_guess.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/gsi/prt_guess.f90 b/src/gsi/prt_guess.f90 index 11befcd2da..469cddbf21 100644 --- a/src/gsi/prt_guess.f90 +++ b/src/gsi/prt_guess.f90 @@ -341,6 +341,7 @@ subroutine prt_guess2(sgrep) ! get pointer to cloud water condensate ier=0;nvarsc=0 iql=0;iqi=0;iqr=0;iqs=0;iqg=0 + ini=0;inr=0 call gsi_metguess_get('clouds::3d',n_actual_clouds,ier) if (mype==0) write(6,*)'prt_guess2: n_actual_clouds = ', n_actual_clouds if (n_actual_clouds>0) then @@ -387,7 +388,7 @@ subroutine prt_guess2(sgrep) if (ier==0) nvarsc=nvarsc+1 endif end if - nvarsc=n_actual_clouds + !nvarsc=n_actual_clouds call gsi_metguess_get ( 'var::cf', ivar, ier ); icf=ivar if (ivar > 0) then call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'cf',ges_cf_it,istatus) From f91f63f89979c2841f7cfbac946c5d2997a04097 Mon Sep 17 00:00:00 2001 From: Azadeh Date: Tue, 21 May 2024 03:34:54 +0000 Subject: [PATCH 06/13] Applied Emily's changes --- src/gsi/crtm_interface.f90 | 42 ++++++++++----------------------- src/gsi/general_read_gfsatm.f90 | 5 ++-- src/gsi/prt_guess.f90 | 2 +- src/gsi/set_crtm_cloudmod.f90 | 2 +- 4 files changed, 18 insertions(+), 33 deletions(-) diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90 index d095c212c9..fd2949b6b0 100644 --- a/src/gsi/crtm_interface.f90 +++ b/src/gsi/crtm_interface.f90 @@ -864,13 +864,9 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo endif ! nvege_type endif ! regional or IGBP -! Initial GFDL saturation water vapor pressure tables - -!Option for Thompson microphysics - +! Initial GFDL saturation water vapor pressure tables if (use_gfdl_qsat) then - if (n_actual_aerosols_wk>0 .or. n_clouds_fwd_wk>0 .and. imp_physics==11 .or. imp_physics==08) then - + if (n_actual_aerosols_wk>0 .or. n_clouds_fwd_wk>0 .and. imp_physics==11) then if (mype==0) write(6,*)myname_,':initial and load GFDL saturation water vapor pressure tables' allocate(table (length)) @@ -1134,7 +1130,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & reshape((/0.0_r_kind, 1.0_r_kind, 1.0_r_kind, 2.0_r_kind, 1.0_r_kind, & -1.0_r_kind, 1.0_r_kind, -1.0_r_kind/), (/4, 2/)) real(r_kind),parameter:: jac_pert = 1.0_r_kind - + real(r_kind),parameter:: xrc3 = 200.0_r_kind ! Declare local variables integer(i_kind):: iquadrant integer(i_kind):: ier,ii,kk,kk2,i,itype,leap_day,day_of_year @@ -1148,9 +1144,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & integer(i_kind):: error_status_clr integer(i_kind):: idx700,dprs,dprs_min integer(i_kind),dimension(8)::obs_time,anal_time - integer(i_kind),dimension(msig) :: klevel - real(r_kind),dimension(nsig) :: qsat - real(r_kind),dimension(nsig) :: qcond + integer(i_kind),dimension(msig) :: klevel ! ****************************** ! Constrained indexing for lai @@ -1183,6 +1177,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & real(r_kind),dimension(nsig) :: rho_air ! density of air (kg/m3) real(r_kind),dimension(nsig) :: cf_calc ! GFDL cloud fraction calculation real(r_kind),dimension(nsig) :: qmix ! water vapor mixing ratio + real(r_kind),dimension(nsig) :: qcond real(r_kind),allocatable,dimension(:,:) :: tgas1d real(r_kind),pointer,dimension(:,: )::psges_itsig =>NULL() real(r_kind),pointer,dimension(:,: )::psges_itsigp=>NULL() @@ -1205,8 +1200,6 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & real(r_kind),pointer,dimension(:,:,:)::nrges_itsig =>NULL() real(r_kind),pointer,dimension(:,:,:)::nrges_itsigp=>NULL() logical :: lmfdeep2 - real(r_kind) , parameter:: xrc3=200.0_r_kind - logical :: sea,icmask integer(i_kind),parameter,dimension(12):: mday=(/0,31,59,90,& @@ -1214,7 +1207,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & real(r_kind) :: lai m1=mype+1 - + if (mype==0) write(6,*) myname_, ' imp_physics = ', imp_physics if (n_clouds_fwd_wk>0) hwp_guess=zero hwp_total=zero theta_700=zero @@ -1374,9 +1367,9 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & iqs=istatus call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'q',qges_itsigp,istatus) iqs=iqs+istatus - call gsi_bundlegetpointer(gsi_metguess_bundle(itsig ),'cf',cfges_itsig ,icfs) + call gsi_bundlegetpointer(gsi_metguess_bundle(itsig ),'cf',cfges_itsig ,istatus) icfs=istatus - call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'cf',cfges_itsigp,icfs) + call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'cf',cfges_itsigp,istatus) icfs=icfs+istatus call gsi_bundlegetpointer(gsi_metguess_bundle(itsig ),'ni',niges_itsig,istatus) inis=istatus @@ -1761,9 +1754,10 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & ges_qsat (ixp,iy ,k, itsigp)*w10+ & ges_qsat (ix ,iyp,k, itsigp)*w01+ & ges_qsat (ixp,iyp,k, itsigp)*w11)*dtsigp - rh(k) = q(k)/qsat(k) !added for cloud fraction computation + c3(k)=r1000/(one-q(k)) qmix(k)=q(k)*c3(k) !convert specific humidity to mixing ratio + rh(k) = q(k)/qs(k) !calculate relative humidity; added for cloud fraction computation ! Space-time interpolation of ozone(poz) if (iozs==0) then poz(k)=((ozges_itsig (ix ,iy ,k)*w00+ & @@ -1790,7 +1784,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & cfges_itsigp(ix ,iyp,k)*w01+ & cfges_itsigp(ixp,iyp,k)*w11)*dtsigp) -! Ensure ozone is greater than ozsmall +! Ensure cloud fraction is between 0 and 1 cf(k)=min(max(zero,cf(k)),one) endif ! cf @@ -1903,7 +1897,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & cf_calc = zero ! sum of the cloud condensate amount for all 5 hydrometeos (ql + qi + qs + qg + qh + qr) qcond(:) = cloud(:,1) + cloud(:,2) + cloud(:,3) + cloud(:,4) + cloud(:,5) - call calc_thompson_cloudfrac (nsig , lmfdeep2, xrc3, prsl, qcond(:), rh, qsat,cf_calc) + call calc_thompson_cloudfrac (nsig , lmfdeep2, xrc3, prsl, qcond(:), rh, qs, cf_calc) cf = cf_calc icfs = 0 ! load cloud fraction into CRTM endif @@ -2122,21 +2116,11 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & do ii=1,n_clouds_fwd_wk !cloud_cont(k,ii)=cloud(kk2,ii)*kgkg_kgm2 cloud_cont(k,ii)=cloud(kk2,ii)*c6(k) - !Option for Thompson added - if (imp_physics==11 .or. imp_physics==08 .and. lprecip_wk .and. cloud_cont(k,ii) > 1.0e-6_r_kind) then + if (lprecip_wk .and. cloud_cont(k,ii) > 1.0e-6_r_kind) then cloud_efr(k,ii)=cloudefr(kk2,ii) else cloud_efr(k,ii)=zero endif - if (trim(cloud_names_fwd(ii))=='ni' .and. atmosphere(1)%temperature(k) Date: Thu, 30 May 2024 21:26:54 +0000 Subject: [PATCH 07/13] Updated the crtm_interface code for thompson reff subroutine --- src/gsi/crtm_interface.f90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90 index fd2949b6b0..98d8e8f63d 100644 --- a/src/gsi/crtm_interface.f90 +++ b/src/gsi/crtm_interface.f90 @@ -2743,7 +2743,7 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) real(r_kind), parameter :: reff_g_min = 150.0_r_kind ! [micron] ! previous value was 0.0_r_kind real(r_kind), parameter :: reff_g_max = 5000.0_r_kind ! [micron] ! previous value was 10000.0_r_kind real(r_kind), parameter:: mu_g = 0.0_r_kind - real(r_kind), parameter :: no_exp = 1.0e-4_r_kind + real(r_kind), parameter :: no_exp = 10.0_r_kind !cloud water if (trim(cloud_name)=='ql') then @@ -2752,7 +2752,7 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) reff_max = reff_w_max do k = 1, nsig qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) - if (qx < qmin) then + if (qx > qmin) then mu_w = MAX(2, MIN((NINT(1000.E6_r_kind/Nt_c) + 2), 15)) lam_w=exp(1.0_r_kind / 3.0_r_kind * log ((am_w*Nt_c *gamma(mu_w + 3.0_r_kind + 1.0_r_kind))/(qx*gamma(mu_w+1.0_r_kind)))) @@ -2770,7 +2770,7 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) reff_max = reff_i_max do k = 1, nsig qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) - if (qx < qmin) then + if (qx > qmin) then lam_i=exp(1.0_r_kind / 3.0_r_kind * log((am_i*ni(k) *gamma(mu_i + 3.0_r_kind + 1.0_r_kind))/(qx*gamma(mu_i+1.0_r_kind)))) reff(k) = 0.5_r_kind * (3.0_r_kind /lam_i)*1.0e6_r_kind reff(k) = max(reff_min, min(reff_max, reff(k))) @@ -2785,7 +2785,7 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) reff_max = reff_r_max do k = 1, nsig qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) - if (qx < qmin) then + if (qx > qmin) then lam_r=exp(1.0_r_kind / 3.0_r_kind * log ((am_r*nr(k) *gamma(mu_r + 3.0_r_kind + 1.0_r_kind))/(qx*gamma(mu_r + 1.0_r_kind)))) reff(k) = 0.5_r_kind *(3.0_r_kind/lam_r)*1.0e6_r_kind reff(k) = max(reff_min, min(reff_max, reff(k))) @@ -2803,7 +2803,7 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) !..Calculate bm_s+1 (th) moment, smoc. Useful for diameter calcs. tc0 = MIN(-0.1_r_kind, tsen(k)-273.15_r_kind) qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) - if (qx < qmin) then + if (qx > qmin) then smob =qx/am_s loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) & + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 & @@ -2830,7 +2830,7 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) am_g = rho_g*pi/6.0_r_kind do k = 1, nsig qx = qxmr(k)*rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) - if (qx < qmin) then + if (qx > qmin) then lam_exp = no_exp* exp(1.0_r_kind/(3.0_r_kind+1.0_r_kind) * log ((am_g*gamma(3.0_r_kind +1.0_r_kind))/qx)) lam_g=lam_exp*exp(1/3.0_r_kind*log(gamma(3.0_r_kind+mu_g+1.0_r_kind)/((3.0_r_kind+mu_g+1.0_r_kind)*(mu_g+1.0_r_kind)))) reff(k) = 0.5_r_kind *(3.0_r_kind/ lam_g)*1.0e6_r_kind From 59b19a856c4a425b4cd639d6766cfa774627797a Mon Sep 17 00:00:00 2001 From: Azadeh Date: Thu, 6 Jun 2024 15:23:00 +0000 Subject: [PATCH 08/13] adding lmfdee2 = .true. for cloud fraction calculation --- src/gsi/crtm_interface.f90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90 index 98d8e8f63d..09caebb7a5 100644 --- a/src/gsi/crtm_interface.f90 +++ b/src/gsi/crtm_interface.f90 @@ -1894,6 +1894,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & ! Calculate Thompson cloud fraction if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==08 .and. lprecip_wk ) then + lmfdeep2 = .true. cf_calc = zero ! sum of the cloud condensate amount for all 5 hydrometeos (ql + qi + qs + qg + qh + qr) qcond(:) = cloud(:,1) + cloud(:,2) + cloud(:,3) + cloud(:,4) + cloud(:,5) @@ -3102,8 +3103,8 @@ subroutine calc_thompson_cloudfrac(nsig, lmfdeep2, xrc3, prsl, clwf, rhly, qstl, real(r_kind) , dimension(nsig), intent(inout) :: cldtot ! --- local variables: - real(r_kind) :: clwmin, clwm, clwt, onemrh, value, tem1, tem2 - integer(i_kind) :: k + real(r_kind) :: clwmin, clwm, clwt, onemrh, value, tem1, tem2 + integer(i_kind) :: k !> - Compute layer cloud fraction. clwmin = 0.0_r_kind From 28bac26fb7576cb0bac075edff84d2bc6537d689 Mon Sep 17 00:00:00 2001 From: Azadeh Date: Wed, 17 Jul 2024 19:43:30 +0000 Subject: [PATCH 09/13] updated the fix submodule hash to point at 2a86d5b --- fix | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fix b/fix index 15ffa60307..2a86d5b1a0 160000 --- a/fix +++ b/fix @@ -1 +1 @@ -Subproject commit 15ffa60307bbc19746d8caeb41782de0b7604be6 +Subproject commit 2a86d5b1a09a8908618d1c8a376c68e8d9b3d02c From 0181f4576a358cca0a0b33c4c51c7a3b7e98ab9f Mon Sep 17 00:00:00 2001 From: Azadeh Date: Thu, 18 Jul 2024 15:39:08 +0000 Subject: [PATCH 10/13] allow precipitation in EnkF update --- src/enkf/gridio_gfs.f90 | 6 ------ 1 file changed, 6 deletions(-) diff --git a/src/enkf/gridio_gfs.f90 b/src/enkf/gridio_gfs.f90 index 35a0c3fbe4..6e68f1a347 100644 --- a/src/enkf/gridio_gfs.f90 +++ b/src/enkf/gridio_gfs.f90 @@ -195,9 +195,6 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & sst_ind = getindex(vars2d, 'sst') use_full_hydro = ( ql_ind > 0 .and. qi_ind > 0 .and. & qr_ind > 0 .and. qs_ind > 0 .and. qg_ind > 0 ) - ! Currently, we do not let precipiation to affect the enkf analysis - ! The following line will be removed after testing - use_full_hydro = .false. if (.not. isinitialized) call init_spec_vars(nlons,nlats,ntrunc,4) @@ -459,9 +456,6 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & end if ! cloud derivatives - ! Currently, we do not let precipiation to affect the enkf analysis - ! The following line will be removed after testing - use_full_hydro = .true. if (.not. use_full_hydro .and. iope==0) then if (ql_ind > 0 .or. qi_ind > 0) then do k=1,nlevs From 3bf46cb6da3477b3f0bdc3eea4979f6d676d4392 Mon Sep 17 00:00:00 2001 From: Azadeh Date: Wed, 24 Jul 2024 16:50:24 +0000 Subject: [PATCH 11/13] followed GSI Code Standards for spaces and indentation --- src/gsi/crtm_interface.f90 | 113 ++++++++++++++++++------------------- src/gsi/prt_guess.f90 | 1 - 2 files changed, 55 insertions(+), 59 deletions(-) diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90 index 09caebb7a5..2914445ac5 100644 --- a/src/gsi/crtm_interface.f90 +++ b/src/gsi/crtm_interface.f90 @@ -866,7 +866,7 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo ! Initial GFDL saturation water vapor pressure tables if (use_gfdl_qsat) then - if (n_actual_aerosols_wk>0 .or. n_clouds_fwd_wk>0 .and. imp_physics==11) then + if (n_actual_aerosols_wk>0 .or. n_clouds_fwd_wk>0 .and. imp_physics==11) then if (mype==0) write(6,*)myname_,':initial and load GFDL saturation water vapor pressure tables' allocate(table (length)) @@ -1144,8 +1144,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & integer(i_kind):: error_status_clr integer(i_kind):: idx700,dprs,dprs_min integer(i_kind),dimension(8)::obs_time,anal_time - integer(i_kind),dimension(msig) :: klevel - + integer(i_kind),dimension(msig) :: klevel ! ****************************** ! Constrained indexing for lai ! CRTM 2.1 implementation change @@ -1897,8 +1896,8 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & lmfdeep2 = .true. cf_calc = zero ! sum of the cloud condensate amount for all 5 hydrometeos (ql + qi + qs + qg + qh + qr) - qcond(:) = cloud(:,1) + cloud(:,2) + cloud(:,3) + cloud(:,4) + cloud(:,5) - call calc_thompson_cloudfrac (nsig , lmfdeep2, xrc3, prsl, qcond(:), rh, qs, cf_calc) + qcond(:) = cloud(:,1) + cloud(:,2) + cloud(:,3) + cloud(:,4) + cloud(:,5) + call calc_thompson_cloudfrac (nsig , lmfdeep2, xrc3, prsl, qcond(:), rh, qs, cf_calc) cf = cf_calc icfs = 0 ! load cloud fraction into CRTM endif @@ -2754,11 +2753,11 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) do k = 1, nsig qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) if (qx > qmin) then - mu_w = MAX(2, MIN((NINT(1000.E6_r_kind/Nt_c) + 2), 15)) - lam_w=exp(1.0_r_kind / 3.0_r_kind * log ((am_w*Nt_c *gamma(mu_w + 3.0_r_kind + 1.0_r_kind))/(qx*gamma(mu_w+1.0_r_kind)))) + mu_w = MAX(2, MIN((NINT(1000.E6_r_kind/Nt_c) + 2), 15)) + lam_w=exp(1.0_r_kind / 3.0_r_kind * log ((am_w*Nt_c *gamma(mu_w + 3.0_r_kind + 1.0_r_kind))/(qx*gamma(mu_w+1.0_r_kind)))) - reff(k) = 0.5_r_kind * ((3.0_r_kind+mu_w)/lam_w)*1.0e6_r_kind - reff(k) = max(reff_min, min(reff_max, reff(k))) + reff(k) = 0.5_r_kind * ((3.0_r_kind+mu_w)/lam_w)*1.0e6_r_kind + reff(k) = max(reff_min, min(reff_max, reff(k))) else reff(k) = zero endif @@ -2772,11 +2771,11 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) do k = 1, nsig qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) if (qx > qmin) then - lam_i=exp(1.0_r_kind / 3.0_r_kind * log((am_i*ni(k) *gamma(mu_i + 3.0_r_kind + 1.0_r_kind))/(qx*gamma(mu_i+1.0_r_kind)))) - reff(k) = 0.5_r_kind * (3.0_r_kind /lam_i)*1.0e6_r_kind - reff(k) = max(reff_min, min(reff_max, reff(k))) + lam_i=exp(1.0_r_kind / 3.0_r_kind * log((am_i*ni(k) *gamma(mu_i + 3.0_r_kind + 1.0_r_kind))/(qx*gamma(mu_i+1.0_r_kind)))) + reff(k) = 0.5_r_kind * (3.0_r_kind /lam_i)*1.0e6_r_kind + reff(k) = max(reff_min, min(reff_max, reff(k))) else - reff(k) = zero + reff(k) = zero endif enddo !Rain @@ -2788,8 +2787,8 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) if (qx > qmin) then lam_r=exp(1.0_r_kind / 3.0_r_kind * log ((am_r*nr(k) *gamma(mu_r + 3.0_r_kind + 1.0_r_kind))/(qx*gamma(mu_r + 1.0_r_kind)))) - reff(k) = 0.5_r_kind *(3.0_r_kind/lam_r)*1.0e6_r_kind - reff(k) = max(reff_min, min(reff_max, reff(k))) + reff(k) = 0.5_r_kind *(3.0_r_kind/lam_r)*1.0e6_r_kind + reff(k) = max(reff_min, min(reff_max, reff(k))) else reff(k) = zero endif @@ -2805,20 +2804,20 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) tc0 = MIN(-0.1_r_kind, tsen(k)-273.15_r_kind) qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) if (qx > qmin) then - smob =qx/am_s - loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) & - + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 & - + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) & - + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 & - + sa(10)*cse(1)*cse(1)*cse(1) - a_ = 10.0**loga_ - b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) & + smob =qx/am_s + loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) & + + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 & + + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) & + + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 & + + sa(10)*cse(1)*cse(1)*cse(1) + a_ = 10.0**loga_ + b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) & + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) & + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1) smoc = a_ * smob**b_ reff(k) = MAX(2.51E-6_r_kind, MIN(0.5*(smoc/smob), 1999.E-6_r_kind)) - reff(k) = max(reff_min, min(reff_max, reff(k)*1.0e6_r_kind)) + reff(k) = max(reff_min, min(reff_max, reff(k)*1.0e6_r_kind)) else reff(k) = zero endif @@ -3094,47 +3093,45 @@ subroutine calc_thompson_cloudfrac(nsig, lmfdeep2, xrc3, prsl, clwf, rhly, qstl, implicit none ! --- inputs: - integer(i_kind), intent(in) :: nsig - real(r_kind) , dimension(nsig), intent(in) :: prsl, clwf, rhly, qstl - logical, intent(in) :: lmfdeep2 - real(r_kind) , intent(in):: xrc3 + integer(i_kind), intent(in) :: nsig + real(r_kind) , dimension(nsig), intent(in) :: prsl, clwf, rhly, qstl + logical, intent(in) :: lmfdeep2 + real(r_kind) , intent(in):: xrc3 ! --- outputs - real(r_kind) , dimension(nsig), intent(inout) :: cldtot + real(r_kind) , dimension(nsig), intent(inout) :: cldtot ! --- local variables: - real(r_kind) :: clwmin, clwm, clwt, onemrh, value, tem1, tem2 - integer(i_kind) :: k - + real(r_kind) :: clwmin, clwm, clwt, onemrh, value, tem1, tem2 + integer(i_kind) :: k + !> - Compute layer cloud fraction. - clwmin = 0.0_r_kind - do k = 1, nsig-1 + clwmin = 0.0_r_kind + do k = 1, nsig-1 ! converting prsl from cb to bar - clwt = 1.0e-10_r_kind * (prsl(k)*0.01_r_kind) - if (clwf(k) > clwt) then - if(rhly(k) > 0.99_r_kind) then - cldtot(k) = 1.0_r_kind - else - onemrh= max( 1.e-10_r_kind, 1.0_r_kind-rhly(k) ) - clwm = clwmin / max( 0.01_r_kind, prsl(k)*0.01_r_kind ) - - tem1 = min(max((onemrh*qstl(k))**0.49_r_kind,0.0001_r_kind),1.0_r_kind) - if (lmfdeep2) then - tem1 = xrc3 / tem1 - else - tem1 = 100.0_r_kind / tem1 - endif - - value = max( min( tem1*(clwf(k)-clwm), 50.0_r_kind ), 0.0_r_kind ) - tem2 = sqrt( sqrt(rhly(k)) ) - cldtot(k) = max( tem2*(1.0_r_kind-exp(-value)), 0.0_r_kind ) - endif - else - cldtot(k) = 0.0_r_kind - endif - enddo + clwt = 1.0e-10_r_kind * (prsl(k)*0.01_r_kind) + if (clwf(k) > clwt) then + if(rhly(k) > 0.99_r_kind) then + cldtot(k) = 1.0_r_kind + else + onemrh= max( 1.e-10_r_kind, 1.0_r_kind-rhly(k) ) + clwm = clwmin / max( 0.01_r_kind, prsl(k)*0.01_r_kind ) + tem1 = min(max((onemrh*qstl(k))**0.49_r_kind,0.0001_r_kind),1.0_r_kind) + if (lmfdeep2) then + tem1 = xrc3 / tem1 + else + tem1 = 100.0_r_kind / tem1 + endif + value = max( min( tem1*(clwf(k)-clwm), 50.0_r_kind ), 0.0_r_kind ) + tem2 = sqrt( sqrt(rhly(k)) ) + cldtot(k) = max( tem2*(1.0_r_kind-exp(-value)), 0.0_r_kind ) + endif + else + cldtot(k) = 0.0_r_kind + endif + enddo - end subroutine calc_thompson_cloudfrac + end subroutine calc_thompson_cloudfrac !............................................ subroutine qs_table(n) diff --git a/src/gsi/prt_guess.f90 b/src/gsi/prt_guess.f90 index 99e5cfd6a2..dbce7e0f79 100644 --- a/src/gsi/prt_guess.f90 +++ b/src/gsi/prt_guess.f90 @@ -388,7 +388,6 @@ subroutine prt_guess2(sgrep) if (ier==0) nvarsc=nvarsc+1 endif end if - !nvarsc=n_actual_clouds call gsi_metguess_get ( 'var::cf', ivar, ier ); icf=ivar if (ivar > 0) then call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'cf',ges_cf_it,istatus) From e28d668041e29884b0da7c2aa08238101845357b Mon Sep 17 00:00:00 2001 From: Azadeh Date: Fri, 26 Jul 2024 16:15:21 +0000 Subject: [PATCH 12/13] removed unused variables --- src/gsi/crtm_interface.f90 | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90 index 2914445ac5..5de662b91d 100644 --- a/src/gsi/crtm_interface.f90 +++ b/src/gsi/crtm_interface.f90 @@ -2694,11 +2694,10 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) ! Declare local variables character(len=*), parameter :: myname_ = 'calc_thompson_reff' integer(i_kind) :: k - integer(i_kind) :: j integer(i_kind) :: mu_w real(r_kind) :: qx real(r_kind) :: reff_min, reff_max - real(r_kind):: lam_i,lam_w,lam_s,lam_r, lam_g,lam_exp, am_r,am_w,am_i,am_g, N_0r ,N_0i,N_0w + real(r_kind):: lam_i,lam_w,lam_r, lam_g,lam_exp, am_r,am_w,am_i,am_g ! Parameters real(r_kind), parameter :: qmin = 1.0e-12_r_kind ! [kg/kg ] @@ -2736,7 +2735,7 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) real(r_kind), parameter :: am_s = 0.069_r_kind real(r_kind), parameter :: bm_s = 2.0_r_kind real(r_kind), dimension(1), parameter :: cse = (/ bm_s + 1.0_r_kind /) - real(r_kind) :: tc0, smo2, smob, smoc, smoz, a_, b_, loga_ + real(r_kind) :: tc0, smob, smoc, a_, b_, loga_ ! Parameters for graupel (Lin 1983) real(r_kind), parameter :: rho_g = 500.0_r_kind ! [kg/m3 ] @@ -2810,7 +2809,7 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) & + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 & + sa(10)*cse(1)*cse(1)*cse(1) - a_ = 10.0**loga_ + a_ = 10.0_r_kind**loga_ b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) & + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) & + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) & From 1daf07bee2d2a18d542fdd31d044d52e3cecb4aa Mon Sep 17 00:00:00 2001 From: Azadeh Date: Fri, 26 Jul 2024 16:56:38 +0000 Subject: [PATCH 13/13] add _r_kind to 0.5 --- src/gsi/crtm_interface.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90 index 5de662b91d..32a2cdb067 100644 --- a/src/gsi/crtm_interface.f90 +++ b/src/gsi/crtm_interface.f90 @@ -2815,7 +2815,7 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1) smoc = a_ * smob**b_ - reff(k) = MAX(2.51E-6_r_kind, MIN(0.5*(smoc/smob), 1999.E-6_r_kind)) + reff(k) = MAX(2.51E-6_r_kind, MIN(0.5_r_kind*(smoc/smob), 1999.E-6_r_kind)) reff(k) = max(reff_min, min(reff_max, reff(k)*1.0e6_r_kind)) else reff(k) = zero