diff --git a/src/enkf/gridio_gfs.f90 b/src/enkf/gridio_gfs.f90 index fe5199e39..456f8126e 100644 --- a/src/enkf/gridio_gfs.f90 +++ b/src/enkf/gridio_gfs.f90 @@ -157,7 +157,7 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & levdim = get_dim(dset,'pfull'); nlevsin = levdim%len idvc=2 else - print *, 'parallel read only supported for netCDF, stopping with error' + write(6,*)'READGRIDDATA_PNC: ***FATAL ERROR*** parallel read only supported for netCDF' , ' PROGRAM STOPS' call stop2(23) end if ice = .false. ! calculate qsat w/resp to ice? @@ -185,6 +185,9 @@ 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) @@ -195,7 +198,7 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & call read_vardata(dset, 'pressfc', values_2d,errcode=iret) if (iret /= 0) then - print *,'error reading ps' + write(6,*)'READGRIDDATA_PNC: ***FATAL ERROR*** reading ps, iret= ',iret,' PROGRAM STOPS' call stop2(31) endif psg = 0.01_r_kind*reshape(values_2d,(/nlons*nlats/)) @@ -215,12 +218,12 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & call read_vardata(dset, 'ugrd', ug3d, ncstart=ncstart, nccount=nccount, errcode=iret) if (iret /= 0) then - print *,'error reading ugrd' + write(6,*)'READGRIDDATA_PNC: ***FATAL ERROR*** reading ugrd, iret= ',iret,' PROGRAM STOPS' call stop2(22) endif call read_vardata(dset, 'vgrd', vg3d, ncstart=ncstart, nccount=nccount, errcode=iret) if (iret /= 0) then - print *,'error reading vgrd' + write(6,*)'READGRIDDATA_PNC: ***FATAL ERROR*** reading vgrd, iret= ',iret,' PROGRAM STOPS' call stop2(23) endif call mpi_gatherv(ug3d, recvcounts(iope+1), mpi_real4, ug3d_0, recvcounts, displs,& @@ -247,12 +250,12 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & end if call read_vardata(dset,'tmp', ug3d, ncstart=ncstart, nccount=nccount, errcode=iret) if (iret /= 0) then - print *,'error reading tmp' + write(6,*)'READGRIDDATA_PNC: ***FATAL ERROR*** reading tmp, iret= ',iret,' PROGRAM STOPS' call stop2(24) endif call read_vardata(dset,'spfh', vg3d, ncstart=ncstart, nccount=nccount, errcode=iret) if (iret /= 0) then - print *,'error reading spfh' + write(6,*)'READGRIDDATA_PNC: ***FATAL ERROR*** reading spfh, iret= ',iret,' PROGRAM STOPS' call stop2(25) endif call mpi_gatherv(ug3d, recvcounts(iope+1), mpi_real4, ug3d_0, recvcounts, displs,& @@ -276,7 +279,7 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & if (oz_ind > 0) then call read_vardata(dset, 'o3mr', ug3d, ncstart=ncstart, nccount=nccount, errcode=iret) if (iret /= 0) then - print *,'error reading o3mr' + write(6,*)'READGRIDDATA_PNC: ***FATAL ERROR*** reading o3mr, iret= ',iret,' PROGRAM STOPS' call stop2(26) endif if (cliptracers) where (ug3d < clip) ug3d = clip @@ -290,31 +293,122 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & end do end if endif - if (cw_ind > 0 .or. ql_ind > 0 .or. qi_ind > 0) then - call read_vardata(dset, 'clwmr', ug3d, ncstart=ncstart, nccount=nccount, errcode=iret) - if (iret /= 0) then - print *,'error reading clwmr' - call stop2(27) + ! Read in hydrometeor fields based on control/state variables listed in anavinfo table + if (use_full_hydro) then + if(ql_ind > 0) then + call read_vardata(dset, 'clwmr', ug3d, ncstart=ncstart, nccount=nccount, errcode=iret) + if (iret /= 0) then + write(6,*)'READGRIDDATA_PNC: ***FATAL ERROR*** reading clwmr, iret= ',iret,' PROGRAM STOPS' + call stop2(26) + endif + if (cliptracers) where (ug3d < clip) ug3d = clip + call mpi_gatherv(ug3d, recvcounts(iope+1), mpi_real4, ug3d_0, recvcounts, displs,& + mpi_real4, 0, iocomms(mem_pe(nproc)),iret) + if (iope==0) then + do k=1,nlevs + krev = nlevs-k+1 + ug = reshape(ug3d_0(:,:,krev),(/nlons*nlats/)) + call copytogrdin(ug,grdin(:,levels(ql_ind-1)+k,nb,ne)) + end do + end if endif - if (imp_physics == 11) then - call read_vardata(dset, 'icmr', vg3d, ncstart=ncstart, nccount=nccount, errcode=iret) + if(qi_ind > 0) then + call read_vardata(dset, 'icmr', ug3d, ncstart=ncstart, nccount=nccount, errcode=iret) if (iret /= 0) then - print *,'error reading icmr' - call stop2(28) + write(6,*)'READGRIDDATA_PNC: ***FATAL ERROR*** reading icmr, iret= ',iret,' PROGRAM STOPS' + call stop2(26) endif - ug3d = ug3d + vg3d + if (cliptracers) where (ug3d < clip) ug3d = clip + call mpi_gatherv(ug3d, recvcounts(iope+1), mpi_real4, ug3d_0, recvcounts, displs,& + mpi_real4, 0, iocomms(mem_pe(nproc)),iret) + if (iope==0) then + do k=1,nlevs + krev = nlevs-k+1 + ug = reshape(ug3d_0(:,:,krev),(/nlons*nlats/)) + call copytogrdin(ug,grdin(:,levels(qi_ind-1)+k,nb,ne)) + end do + end if + endif + if(qr_ind > 0) then + call read_vardata(dset, 'rwmr', ug3d, ncstart=ncstart, nccount=nccount, errcode=iret) + if (iret /= 0) then + write(6,*)'READGRIDDATA_PNC: ***FATAL ERROR*** reading rwmr, iret= ',iret,' PROGRAM STOPS' + call stop2(26) + endif + if (cliptracers) where (ug3d < clip) ug3d = clip + call mpi_gatherv(ug3d, recvcounts(iope+1), mpi_real4, ug3d_0, recvcounts, displs,& + mpi_real4, 0, iocomms(mem_pe(nproc)),iret) + if (iope==0) then + do k=1,nlevs + krev = nlevs-k+1 + ug = reshape(ug3d_0(:,:,krev),(/nlons*nlats/)) + call copytogrdin(ug,grdin(:,levels(qr_ind-1)+k,nb,ne)) + end do + end if + endif + if(qs_ind > 0) then + call read_vardata(dset, 'snmr', ug3d, ncstart=ncstart, nccount=nccount, errcode=iret) + if (iret /= 0) then + write(6,*)'READGRIDDATA_PNC: ***FATAL ERROR*** reading snmr, iret= ',iret,' PROGRAM STOPS' + call stop2(26) + endif + if (cliptracers) where (ug3d < clip) ug3d = clip + call mpi_gatherv(ug3d, recvcounts(iope+1), mpi_real4, ug3d_0, recvcounts, displs,& + mpi_real4, 0, iocomms(mem_pe(nproc)),iret) + if (iope==0) then + do k=1,nlevs + krev = nlevs-k+1 + ug = reshape(ug3d_0(:,:,krev),(/nlons*nlats/)) + call copytogrdin(ug,grdin(:,levels(qs_ind-1)+k,nb,ne)) + end do + end if + endif + if(qg_ind > 0) then + call read_vardata(dset, 'grle', ug3d, ncstart=ncstart, nccount=nccount, errcode=iret) + if (iret /= 0) then + write(6,*)'READGRIDDATA_PNC: ***FATAL ERROR*** reading grle, iret= ',iret,' PROGRAM STOPS' + call stop2(26) + endif + if (cliptracers) where (ug3d < clip) ug3d = clip + call mpi_gatherv(ug3d, recvcounts(iope+1), mpi_real4, ug3d_0, recvcounts, displs,& + mpi_real4, 0, iocomms(mem_pe(nproc)),iret) + if (iope==0) then + do k=1,nlevs + krev = nlevs-k+1 + ug = reshape(ug3d_0(:,:,krev),(/nlons*nlats/)) + call copytogrdin(ug,grdin(:,levels(qg_ind-1)+k,nb,ne)) + end do + end if + endif + else + ! Handle non-precipiting hydrometeors + ! if control or state variable is cw, make sure combine background ql and qi to cw + if (cw_ind > 0 .or. ql_ind > 0 .or. qi_ind > 0) then + call read_vardata(dset, 'clwmr', ug3d, ncstart=ncstart, nccount=nccount, errcode=iret) + if (iret /= 0) then + write(6,*)'READGRIDDATA_PNC: ***FATAL ERROR*** reading clwmr, iret= ',iret,' PROGRAM STOPS' + call stop2(27) + endif + if (imp_physics == 11) then + call read_vardata(dset, 'icmr', vg3d, ncstart=ncstart, nccount=nccount, errcode=iret) + if (iret /= 0) then + write(6,*)'READGRIDDATA_PNC: ***FATAL ERROR*** reading icmr, iret= ',iret,' PROGRAM STOPS' + call stop2(28) + endif + ug3d = ug3d + vg3d + endif + if (cliptracers) where (ug3d < clip) ug3d = clip + call mpi_gatherv(ug3d, recvcounts(iope+1), mpi_real4, ug3d_0, recvcounts, displs,& + mpi_real4, 0, iocomms(mem_pe(nproc)),iret) + if (iope==0) then + do k=1,nlevs + krev = nlevs-k+1 + ug = reshape(ug3d_0(:,:,krev),(/nlons*nlats/)) + call copytogrdin(ug,cw(:,k)) + if (cw_ind > 0) grdin(:,levels(cw_ind-1)+k,nb,ne) = cw(:,k) + end do + end if endif - if (cliptracers) where (ug3d < clip) ug3d = clip - call mpi_gatherv(ug3d, recvcounts(iope+1), mpi_real4, ug3d_0, recvcounts, displs,& - mpi_real4, 0, iocomms(mem_pe(nproc)),iret) - if (iope==0) then - do k=1,nlevs - krev = nlevs-k+1 - ug = reshape(ug3d_0(:,:,krev),(/nlons*nlats/)) - call copytogrdin(ug,cw(:,k)) - if (cw_ind > 0) grdin(:,levels(cw_ind-1)+k,nb,ne) = cw(:,k) - end do - end if endif deallocate(ug3d,vg3d) @@ -355,22 +449,25 @@ 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 - do i = 1, npts - qi_coef = -r0_05*(tv(i,k)/(one+fv*q(i,k))-t0c) - qi_coef = max(zero,qi_coef) - qi_coef = min(one,qi_coef) ! 0<=qi_coef<=1 - if (ql_ind > 0) then - grdin(i,levels(ql_ind-1)+k,nb,ne) = cw(i,k)*(one-qi_coef) - endif - if (qi_ind > 0) then - grdin(i,levels(qi_ind-1)+k,nb,ne) = cw(i,k)*qi_coef - endif + if (ql_ind > 0 .or. qi_ind > 0) then + do k=1,nlevs + do i = 1, npts + qi_coef = -r0_05*(tv(i,k)/(one+fv*q(i,k))-t0c) + qi_coef = max(zero,qi_coef) + qi_coef = min(one,qi_coef) ! 0<=qi_coef<=1 + if (ql_ind > 0) then + grdin(i,levels(ql_ind-1)+k,nb,ne) = cw(i,k)*(one-qi_coef) + endif + if (qi_ind > 0) then + grdin(i,levels(qi_ind-1)+k,nb,ne) = cw(i,k)*qi_coef + endif + enddo enddo - enddo - endif + endif endif if (sst_ind > 0 .and. iope==0) then @@ -3709,7 +3806,7 @@ subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate use netcdf use params, only: nbackgrounds,incfileprefixes,fgfileprefixes,reducedgrid,& datestring,nhr_anal - use constants, only: grav + use constants, only: grav,qcmin use mpi use module_ncio, only: Dataset, Variable, Dimension, open_dataset,& read_attribute, close_dataset, get_dim, read_vardata,& @@ -3744,7 +3841,8 @@ subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate integer(i_kind) :: ncid_out, lon_dimid, lat_dimid, lev_dimid, ilev_dimid integer(i_kind) :: lonvarid, latvarid, levvarid, pfullvarid, ilevvarid, & hyaivarid, hybivarid, uvarid, vvarid, delpvarid, delzvarid, & - tvarid, sphumvarid, liqwatvarid, o3varid, icvarid + tvarid, sphumvarid, liqwatvarid, o3varid, icvarid, & + rwmrvarid, snmrvarid, grlevarid integer(i_kind) :: iadateout ! fixed fields such as lat, lon, levs @@ -3846,6 +3944,12 @@ subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate call nccheck_incr(nf90_var_par_access(ncid_out, o3varid, nf90_collective)) call nccheck_incr(nf90_def_var(ncid_out, "icmr_inc", nf90_real, dimids3, icvarid)) call nccheck_incr(nf90_var_par_access(ncid_out, icvarid, nf90_collective)) + call nccheck_incr(nf90_def_var(ncid_out, "rwmr_inc", nf90_real, dimids3, rwmrvarid)) + call nccheck_incr(nf90_var_par_access(ncid_out, rwmrvarid, nf90_collective)) + call nccheck_incr(nf90_def_var(ncid_out, "snmr_inc", nf90_real, dimids3, snmrvarid)) + call nccheck_incr(nf90_var_par_access(ncid_out, snmrvarid, nf90_collective)) + call nccheck_incr(nf90_def_var(ncid_out, "grle_inc", nf90_real, dimids3, grlevarid)) + call nccheck_incr(nf90_var_par_access(ncid_out, grlevarid, nf90_collective)) ! place global attributes to parallel calc_increment output call nccheck_incr(nf90_put_att(ncid_out, nf90_global, "source", "GSI EnKF")) call nccheck_incr(nf90_put_att(ncid_out, nf90_global, "comment", & @@ -3874,11 +3978,14 @@ subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate ! old logical massbal_adjust, if non-zero 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. dsfg = open_dataset(filenamein, paropen=.true., mpicomm=iocomms(mem_pe(nproc))) call read_attribute(dsfg, 'ak', values_1d,errcode=iret) if (iret /= 0) then - print *,'error reading ak' + write(6,*)'WRITEINCREMENT_PNC: ***FATAL ERROR*** reading ak, iret= ',iret,' PROGRAM STOPS' call stop2(29) endif do k=1,nlevs+1 @@ -3887,7 +3994,7 @@ subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate enddo call read_attribute(dsfg, 'bk', values_1d,errcode=iret) if (iret /= 0) then - print *,'error reading bk' + write(6,*)'WRITEINCREMENT_PNC: ***FATAL ERROR*** reading bk, iret= ',iret,' PROGRAM STOPS' call stop2(29) endif do k=1,nlevs+1 @@ -3997,7 +4104,7 @@ subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate allocate(tvanl(nlons,nlats,nccount(3)),tmpanl(nlons,nlats,nccount(3)),qanl(nlons,nlats,nccount(3))) call read_vardata(dsfg, 'spfh', q, ncstart=ncstart, nccount=nccount, errcode=iret) if (iret /= 0) then - print *,'error reading spfh' + write(6,*)'WRITEINCREMENT_PNC: ***FATAL ERROR*** reading spfh, iret= ',iret,' PROGRAM STOPS' call stop2(29) endif do k=lev_pe1(iope), lev_pe2(iope) @@ -4022,7 +4129,7 @@ subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate ! t increment call read_vardata(dsfg, 'tmp', tmp, ncstart=ncstart, nccount=nccount, errcode=iret) if (iret /= 0) then - print *,'error reading tmp' + write(6,*)'WRITEINCREMENT_PNC: ***FATAL ERROR*** reading tmp, iret= ',iret,' PROGRAM STOPS' call stop2(29) endif tv = tmp * ( 1.0 + fv*q) @@ -4094,39 +4201,171 @@ subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate call nccheck_incr(nf90_put_var(ncid_out, o3varid, sngl(inc3dout), & start = ncstart, count = nccount)) - ! liq wat increment - ! icmr increment - do k=lev_pe1(iope), lev_pe2(iope) - krev = nlevs-k+1 - ki = k - lev_pe1(iope) + 1 - ug = zero - if (cw_ind > 0) then - call copyfromgrdin(grdin(:,levels(cw_ind-1)+krev,nb,ne),ug) - end if - if (imp_physics == 11) then - work = -r0_05 * (reshape(tmpanl(:,:,ki),(/nlons*nlats/)) - t0c) - do i=1,nlons*nlats - work(i) = max(zero,work(i)) - work(i) = min(one,work(i)) - enddo - vg = ug * work ! cloud ice - ug = ug * (one - work) ! cloud water - inc3d2(:,:,ki) = reshape(vg,(/nlons,nlats/)) + ! For hydrometeors, following the treatment for specific humidity increment + ! Need to make sure the analysis value is not negative + ! Read in background + increment and make sure the minimum is qcmin + ! Adjust increment accordingly + + if (use_full_hydro) then + ! liq wat increment + call read_vardata(dsfg, 'clwmr', q, ncstart=ncstart, nccount=nccount, errcode=iret) + if (iret /= 0) then + write(6,*)'WRITEINCREMENT_PNC: ***FATAL ERROR*** reading clwmr, iret= ',iret,' PROGRAM STOPS' + call stop2(29) endif - inc3d(:,:,ki) = reshape(ug,(/nlons,nlats/)) - enddo - do j=1,nlats - inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) - end do - if (should_zero_increments_for('liq_wat_inc')) inc3dout = zero - call nccheck_incr(nf90_put_var(ncid_out, liqwatvarid, sngl(inc3dout), & - start = ncstart, count = nccount)) - do j=1,nlats - inc3dout(:,nlats-j+1,:) = inc3d2(:,j,:) - end do - if (should_zero_increments_for('icmr_inc')) inc3dout = zero - call nccheck_incr(nf90_put_var(ncid_out, icvarid, sngl(inc3dout), & - start = ncstart, count = nccount)) + do k=lev_pe1(iope), lev_pe2(iope) + krev = nlevs-k+1 + ki = k - lev_pe1(iope) + 1 + inc(:) = zero + if (ql_ind > 0) then + call copyfromgrdin(grdin(:,levels(ql_ind-1) + krev,nb,ne),inc) + endif + inc3d(:,:,ki) = reshape(inc,(/nlons,nlats/)) + qanl(:,:,ki) = q(:,:,ki) + inc3d(:,:,ki) + end do + if (cliptracers) where (qanl < qcmin) qanl = qcmin + inc3d = qanl - q ! updated clwmr increment + do j=1,nlats + inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) + end do + if (should_zero_increments_for('liq_wat_inc')) inc3dout = zero + call nccheck_incr(nf90_put_var(ncid_out, liqwatvarid, sngl(inc3dout), & + start = ncstart, count = nccount)) + + ! icmr increment + call read_vardata(dsfg, 'icmr', q, ncstart=ncstart, nccount=nccount, errcode=iret) + if (iret /= 0) then + write(6,*)'WRITEINCREMENT_PNC: ***FATAL ERROR*** reading icmr, iret= ',iret,' PROGRAM STOPS' + call stop2(29) + endif + do k=lev_pe1(iope), lev_pe2(iope) + krev = nlevs-k+1 + ki = k - lev_pe1(iope) + 1 + inc(:) = zero + if (qi_ind > 0) then + call copyfromgrdin(grdin(:,levels(qi_ind-1) + krev,nb,ne),inc) + endif + inc3d(:,:,ki) = reshape(inc,(/nlons,nlats/)) + qanl(:,:,ki) = q(:,:,ki) + inc3d(:,:,ki) + end do + if (cliptracers) where (qanl < qcmin) qanl = qcmin + inc3d = qanl - q ! updated icmr increment + do j=1,nlats + inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) + end do + if (should_zero_increments_for('icmr_inc')) inc3dout = zero + call nccheck_incr(nf90_put_var(ncid_out, icvarid, sngl(inc3dout), & + start = ncstart, count = nccount)) + + ! rwmr increment + call read_vardata(dsfg, 'rwmr', q, ncstart=ncstart, nccount=nccount, errcode=iret) + if (iret /= 0) then + write(6,*)'WRITEINCREMENT_PNC: ***FATAL ERROR*** reading rwmr, iret= ',iret,' PROGRAM STOPS' + call stop2(29) + endif + do k=lev_pe1(iope), lev_pe2(iope) + krev = nlevs-k+1 + ki = k - lev_pe1(iope) + 1 + inc(:) = zero + if (qr_ind > 0) then + call copyfromgrdin(grdin(:,levels(qr_ind-1) + krev,nb,ne),inc) + endif + inc3d(:,:,ki) = reshape(inc,(/nlons,nlats/)) + qanl(:,:,ki) = q(:,:,ki) + inc3d(:,:,ki) + end do + if (cliptracers) where (qanl < qcmin) qanl = qcmin + inc3d = qanl - q ! updated rwmr increment + do j=1,nlats + inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) + end do + if (should_zero_increments_for('rwmr_inc')) inc3dout = zero + call nccheck_incr(nf90_put_var(ncid_out, rwmrvarid, sngl(inc3dout), & + start = ncstart, count = nccount)) + + ! snmr increment + call read_vardata(dsfg, 'snmr', q, ncstart=ncstart, nccount=nccount, errcode=iret) + if (iret /= 0) then + write(6,*)'WRITEINCREMENT_PNC: ***FATAL ERROR*** reading snmr, iret= ',iret,' PROGRAM STOPS' + call stop2(29) + endif + do k=lev_pe1(iope), lev_pe2(iope) + krev = nlevs-k+1 + ki = k - lev_pe1(iope) + 1 + inc(:) = zero + if (qs_ind > 0) then + call copyfromgrdin(grdin(:,levels(qs_ind-1) + krev,nb,ne),inc) + endif + inc3d(:,:,ki) = reshape(inc,(/nlons,nlats/)) + qanl(:,:,ki) = q(:,:,ki) + inc3d(:,:,ki) + end do + if (cliptracers) where (qanl < qcmin) qanl = qcmin + inc3d = qanl - q ! updated snmr increment + do j=1,nlats + inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) + end do + if (should_zero_increments_for('snmr_inc')) inc3dout = zero + call nccheck_incr(nf90_put_var(ncid_out, snmrvarid, sngl(inc3dout), & + start = ncstart, count = nccount)) + + ! grle increment + call read_vardata(dsfg, 'grle', q, ncstart=ncstart, nccount=nccount, errcode=iret) + if (iret /= 0) then + write(6,*)'WRITEINCREMENT_PNC: ***FATAL ERROR*** reading grle, iret= ',iret,' PROGRAM STOPS' + call stop2(29) + endif + do k=lev_pe1(iope), lev_pe2(iope) + krev = nlevs-k+1 + ki = k - lev_pe1(iope) + 1 + inc(:) = zero + if (qg_ind > 0) then + call copyfromgrdin(grdin(:,levels(qg_ind-1) + krev,nb,ne),inc) + endif + inc3d(:,:,ki) = reshape(inc,(/nlons,nlats/)) + qanl(:,:,ki) = q(:,:,ki) + inc3d(:,:,ki) + end do + if (cliptracers) where (qanl < qcmin) qanl = qcmin + inc3d = qanl - q ! updated grle increment + do j=1,nlats + inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) + end do + if (should_zero_increments_for('grle_inc')) inc3dout = zero + call nccheck_incr(nf90_put_var(ncid_out, grlevarid, sngl(inc3dout), & + start = ncstart, count = nccount)) + else + ! liq wat increment + ! icmr increment + do k=lev_pe1(iope), lev_pe2(iope) + krev = nlevs-k+1 + ki = k - lev_pe1(iope) + 1 + ug = zero + if (cw_ind > 0) then + call copyfromgrdin(grdin(:,levels(cw_ind-1)+krev,nb,ne),ug) + end if + if (imp_physics == 11) then + work = -r0_05 * (reshape(tmpanl(:,:,ki),(/nlons*nlats/)) - t0c) + do i=1,nlons*nlats + work(i) = max(zero,work(i)) + work(i) = min(one,work(i)) + enddo + vg = ug * work ! cloud ice + ug = ug * (one - work) ! cloud water + inc3d2(:,:,ki) = reshape(vg,(/nlons,nlats/)) + endif + inc3d(:,:,ki) = reshape(ug,(/nlons,nlats/)) + enddo + do j=1,nlats + inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) + end do + if (should_zero_increments_for('liq_wat_inc')) inc3dout = zero + call nccheck_incr(nf90_put_var(ncid_out, liqwatvarid, sngl(inc3dout), & + start = ncstart, count = nccount)) + do j=1,nlats + inc3dout(:,nlats-j+1,:) = inc3d2(:,j,:) + end do + if (should_zero_increments_for('icmr_inc')) inc3dout = zero + call nccheck_incr(nf90_put_var(ncid_out, icvarid, sngl(inc3dout), & + start = ncstart, count = nccount)) + endif call mpi_barrier(iocomms(mem_pe(nproc)), iret)