diff --git a/src/gsibec/gsi/control2state.F90 b/src/gsibec/gsi/control2state.F90 index 8c1b7a4..a1a3b9b 100644 --- a/src/gsibec/gsi/control2state.F90 +++ b/src/gsibec/gsi/control2state.F90 @@ -105,6 +105,7 @@ subroutine control2state(xhat,sval,bval) integer(i_kind) :: ictd2m,icmxtm,icmitm,icpmsl,ichowv integer(i_kind) :: icsfwter,icvpwter,ictcamt,iclcbas integer(i_kind) :: iccldch,icuwnd10m,icvwnd10m +integer(i_kind) :: icext1,icext2 character(len=3), parameter :: mycvars(ncvars) = (/ & ! vars from CV needed here 'sf ', 'vp ', 'ps ', 't ', 'q ', 'cw ', 'ql ', 'qi ', 'w ' /) logical :: lc_sf,lc_vp,lc_w,lc_ps,lc_t,lc_rh,lc_cw,lc_ql,lc_qi @@ -136,6 +137,7 @@ subroutine control2state(xhat,sval,bval) real(r_kind),pointer,dimension(:,:,:) :: sv_w=>NULL(),sv_dw=>NULL(),sv_prse=>NULL() real(r_kind),pointer,dimension(:,:,:) :: sv_q=>NULL(),sv_tsen=>NULL() real(r_kind),pointer,dimension(:,:,:) :: sv_tv=>NULL(),sv_oz=>NULL() +real(r_kind),pointer,dimension(:,:,:) :: sv_ext1=>NULL(),sv_ext2=>NULL() real(r_kind),pointer,dimension(:,:,:) :: sv_rank3=>NULL() real(r_kind),pointer,dimension(:,:) :: sv_rank2=>NULL() @@ -218,6 +220,8 @@ subroutine control2state(xhat,sval,bval) call gsi_bundlegetpointer (xhat%step(1),'cldch',iccldch,istatus) call gsi_bundlegetpointer (xhat%step(1),'uwnd10m',icuwnd10m,istatus) call gsi_bundlegetpointer (xhat%step(1),'vwnd10m',icvwnd10m,istatus) +call gsi_bundlegetpointer (xhat%step(1),'ext1',icext1,istatus) +call gsi_bundlegetpointer (xhat%step(1),'ext2',icext2,istatus) ! Loop over control steps do jj=1,nsubwin @@ -294,13 +298,13 @@ subroutine control2state(xhat,sval,bval) call gsi_bundlegetpointer (wbundle,'q' ,cv_rh ,istatus) ! Get 3d pressure - if(do_getprs_tl) call getprs_tl(cv_ps,cv_t,sv_prse) +!_RT if(do_getprs_tl) call getprs_tl(cv_ps,cv_t,sv_prse) ! Convert input normalized RH to q if(do_normal_rh_to_q) call normal_rh_to_q(cv_rh,cv_t,sv_prse,sv_q) ! Calculate sensible temperature - if(do_tv_to_tsen .and. .not.regional) call tv_to_tsen(cv_t,sv_q,sv_tsen) +!_RT if(do_tv_to_tsen .and. .not.regional) call tv_to_tsen(cv_t,sv_q,sv_tsen) ! Copy other variables call gsi_bundlegetvar ( wbundle, 't' , sv_tv, istatus ) @@ -411,6 +415,15 @@ subroutine control2state(xhat,sval,bval) call gsi_bundlegetpointer (sval(jj),'vwnd10m' ,sv_vwnd10m, istatus) call gsi_bundlegetvar ( wbundle, 'vwnd10m', sv_vwnd10m, istatus ) end if + if (icext1>0) then + print *, 'DEBUG_RT: cv(tl) ext1' + call gsi_bundlegetpointer (sval(jj),'ext1' ,sv_ext1, istatus) + call gsi_bundlegetvar ( wbundle, 'ext1', sv_ext1, istatus ) + end if + if (icext2>0) then + call gsi_bundlegetpointer (sval(jj),'ext2' ,sv_ext2, istatus) + call gsi_bundlegetvar ( wbundle, 'ext2', sv_ext2, istatus ) + end if ! Same one-to-one map for chemistry-vars; take care of them together do ic=1,ngases diff --git a/src/gsibec/gsi/control2state_ad.F90 b/src/gsibec/gsi/control2state_ad.F90 index 8ce6e71..d5b895c 100644 --- a/src/gsibec/gsi/control2state_ad.F90 +++ b/src/gsibec/gsi/control2state_ad.F90 @@ -97,6 +97,7 @@ subroutine control2state_ad(rval,bval,grad) integer(i_kind) :: ictd2m,icmxtm,icmitm,icpmsl,ichowv integer(i_kind) :: ictcamt,iclcbas,icsfwter,icvpwter integer(i_kind) :: iccldch,icuwnd10m,icvwnd10m +integer(i_kind) :: icext1,icext2 character(len=3), parameter :: mycvars(ncvars) = (/ & 'sf ', 'vp ', 'ps ', 't ', 'q ', 'cw ', 'ql ', 'qi ', 'w ' /) logical :: lc_sf,lc_vp,lc_w,lc_ps,lc_t,lc_rh,lc_cw,lc_ql,lc_qi @@ -126,6 +127,7 @@ subroutine control2state_ad(rval,bval,grad) real(r_kind),pointer,dimension(:,:) :: rv_uwnd10m=>NULL(),rv_vwnd10m=>NULL() real(r_kind),pointer,dimension(:,:,:) :: rv_u=>NULL(),rv_v=>NULL(),rv_w=>NULL(),rv_dw=>NULL(),rv_prse=>NULL() real(r_kind),pointer,dimension(:,:,:) :: rv_q=>NULL(),rv_tsen=>NULL(),rv_tv=>NULL(),rv_oz=>NULL() +real(r_kind),pointer,dimension(:,:,:) :: rv_ext1=>NULL(),rv_ext2=>NULL() real(r_kind),pointer,dimension(:,:,:) :: rv_rank3=>NULL() real(r_kind),pointer,dimension(:,:) :: rv_rank2=>NULL() @@ -204,6 +206,8 @@ subroutine control2state_ad(rval,bval,grad) call gsi_bundlegetpointer (grad%step(1),'cldch',iccldch,istatus) call gsi_bundlegetpointer (grad%step(1),'uwnd10m',icuwnd10m,istatus) call gsi_bundlegetpointer (grad%step(1),'vwnd10m',icvwnd10m,istatus) +call gsi_bundlegetpointer (grad%step(1),'ext1',icext1,istatus) +call gsi_bundlegetpointer (grad%step(1),'ext2',icext2,istatus) ! Loop over control steps do jj=1,nsubwin @@ -306,14 +310,14 @@ subroutine control2state_ad(rval,bval,grad) enddo end if ! Calculate sensible temperature - if(do_tv_to_tsen_ad .and. .not.regional) call tv_to_tsen_ad(cv_t,rv_q,rv_tsen) +!_RT if(do_tv_to_tsen_ad .and. .not.regional) call tv_to_tsen_ad(cv_t,rv_q,rv_tsen) ! Adjoint of convert input normalized RH to q to add contribution of moisture ! to t, p , and normalized rh if(do_normal_rh_to_q_ad) call normal_rh_to_q_ad(cv_rh,cv_t,rv_prse,rv_q) ! Adjoint to convert ps to 3-d pressure - if(do_getprs_ad) call getprs_ad(cv_ps,cv_t,rv_prse) +!_RT if(do_getprs_ad) call getprs_ad(cv_ps,cv_t,rv_prse) !$omp section @@ -413,6 +417,14 @@ subroutine control2state_ad(rval,bval,grad) call gsi_bundlegetpointer (rval(jj),'vwnd10m' ,rv_vwnd10m, istatus) call gsi_bundleputvar ( wbundle, 'vwnd10m', rv_vwnd10m, istatus ) end if + if (icext1>0) then + call gsi_bundlegetpointer (rval(jj),'ext1' ,rv_ext1, istatus) + call gsi_bundleputvar ( wbundle, 'ext1', rv_ext1, istatus ) + end if + if (icext2>0) then + call gsi_bundlegetpointer (rval(jj),'ext2' ,rv_ext2, istatus) + call gsi_bundleputvar ( wbundle, 'ext2', rv_ext2, istatus ) + end if !$omp end parallel sections diff --git a/src/gsibec/gsi/control_vectors.f90 b/src/gsibec/gsi/control_vectors.f90 index b070986..628bf4c 100644 --- a/src/gsibec/gsi/control_vectors.f90 +++ b/src/gsibec/gsi/control_vectors.f90 @@ -118,6 +118,7 @@ module control_vectors ! Public variables ! public cvars2d, cvars3d, cvarsmd, evars2d, evars3d, nrf_var +public fvars2d, fvars3d public nc2d ! number of 2d static control fields public nc3d ! number of 3d static control fields public mvars ! number of motley fields @@ -167,6 +168,8 @@ module control_vectors character(len=max_varname_length),allocatable,dimension(:) :: cvars3d ! 3-d fields for static CV character(len=max_varname_length),allocatable,dimension(:) :: evars2d ! 2-d fields for ensemble CV character(len=max_varname_length),allocatable,dimension(:) :: evars3d ! 3-d fields for ensemble CV +character(len=max_varname_length),allocatable,dimension(:) :: fvars2d ! 2-d CV field names in file +character(len=max_varname_length),allocatable,dimension(:) :: fvars3d ! 3-d CV field names in file character(len=max_varname_length),allocatable,dimension(:) :: cvarsmd ! motley variable names real(r_kind) ,allocatable,dimension(:) :: as3d real(r_kind) ,allocatable,dimension(:) :: as2d @@ -350,6 +353,7 @@ subroutine init_anacv(rcname) enddo allocate(nrf_var(nvars),cvars3d(nc3d),cvars2d(nc2d)) +allocate(fvars3d(nc3d),fvars2d(nc2d)) allocate(as3d(nc3d),as2d(nc2d)) allocate(be3d(nc3d),be2d(nc2d),bemo(mvars)) allocate(cvarsmd(mvars)) @@ -378,11 +382,13 @@ subroutine init_anacv(rcname) if(ilev==1) then nc2d=nc2d+1 cvars2d(nc2d)=trim(adjustl(var)) + fvars2d(nc2d)=trim(adjustl(funcof)) as2d(nc2d)=aas be2d(nc2d)=bes else nc3d=nc3d+1 cvars3d(nc3d)=trim(adjustl(var)) + fvars3d(nc3d)=trim(adjustl(funcof)) nrf_3d(ii)=.true. as3d(nc3d)=aas be3d(nc3d)=bes @@ -422,6 +428,7 @@ subroutine final_anacv deallocate(cvarsmd) deallocate(be3d,be2d,bemo) deallocate(as3d,as2d) + deallocate(fvars2d,fvars3d) deallocate(nrf_var,cvars2d,cvars3d) llinit=.false. end subroutine final_anacv diff --git a/src/gsibec/gsi/ensctl2state.F90 b/src/gsibec/gsi/ensctl2state.F90 index ee651ed..ef0f0f4 100644 --- a/src/gsibec/gsi/ensctl2state.F90 +++ b/src/gsibec/gsi/ensctl2state.F90 @@ -95,6 +95,8 @@ subroutine ensctl2state(xhat,mval,eval) real(r_kind),pointer,dimension(:,:,:) :: sv_rank3=>NULL() real(r_kind),pointer,dimension(:,:,:) :: sv_w=>NULL() real(r_kind),pointer,dimension(:,:,:) :: sv_dw=>NULL() +real(r_kind),pointer,dimension(:,:,:) :: sv_ext1=>NULL() +real(r_kind),pointer,dimension(:,:,:) :: sv_ext2=>NULL() logical :: do_getprs_tl,do_normal_rh_to_q,do_tv_to_tsen,do_getuv,lstrong_bk_vars logical :: do_tlnmc,do_q_copy @@ -185,6 +187,8 @@ subroutine ensctl2state(xhat,mval,eval) call gsi_bundlegetpointer (eval(jj),'u' ,sv_u, istatus) call gsi_bundlegetpointer (eval(jj),'v' ,sv_v, istatus) call gsi_bundlegetpointer (eval(jj),'tsen',sv_tsen,istatus) + call gsi_bundlegetpointer (eval(jj),'ext1',sv_ext1,istatus) + call gsi_bundlegetpointer (eval(jj),'ext2',sv_ext2,istatus) !$omp parallel sections private(ic,id,istatus) !$omp section @@ -207,6 +211,10 @@ subroutine ensctl2state(xhat,mval,eval) ! Copy variables call gsi_bundlegetvar ( wbundle_c, 't' , sv_tv, istatus ) call gsi_bundlegetvar ( wbundle_c, 'ps' , sv_ps, istatus ) + if (associated(sv_ext1)) & + call gsi_bundlegetvar ( wbundle_c, 'ext1',sv_ext1,istatus ) + if (associated(sv_ext2)) & + call gsi_bundlegetvar ( wbundle_c, 'ext2',sv_ext2,istatus ) ! Get 3d pressure if(do_q_copy) then call gsi_bundlegetvar ( wbundle_c, 'q', sv_q, istatus ) diff --git a/src/gsibec/gsi/ensctl2state_ad.F90 b/src/gsibec/gsi/ensctl2state_ad.F90 index ff409a2..d2f4d08 100644 --- a/src/gsibec/gsi/ensctl2state_ad.F90 +++ b/src/gsibec/gsi/ensctl2state_ad.F90 @@ -91,6 +91,8 @@ subroutine ensctl2state_ad(eval,mval,grad) real(r_kind),pointer,dimension(:,:,:) :: rv_rank3=>NULL() real(r_kind),pointer,dimension(:,:,:) :: rv_w=>NULL() real(r_kind),pointer,dimension(:,:,:) :: rv_dw=>NULL() +real(r_kind),pointer,dimension(:,:,:) :: rv_ext1=>NULL() +real(r_kind),pointer,dimension(:,:,:) :: rv_ext2=>NULL() logical :: do_getuv,do_tv_to_tsen_ad,do_normal_rh_to_q_ad,do_getprs_ad,lstrong_bk_vars logical :: do_tlnmc,do_q_copy @@ -174,6 +176,8 @@ subroutine ensctl2state_ad(eval,mval,grad) call gsi_bundlegetpointer (eval(jj),'tsen',rv_tsen,istatus) call gsi_bundlegetpointer (eval(jj),'q' ,rv_q , istatus) call gsi_bundlegetpointer (wbundle_c,'q' ,cv_rh ,istatus) + call gsi_bundlegetpointer (eval(jj),'ext1',rv_ext1 ,istatus) + call gsi_bundlegetpointer (eval(jj),'ext2',rv_ext2 ,istatus) ! Adjoint of consistency for sensible temperature, calculate sensible temperature if(do_tv_to_tsen_ad) call tv_to_tsen_ad(rv_tv,rv_q,rv_tsen) @@ -265,6 +269,10 @@ subroutine ensctl2state_ad(eval,mval,grad) ! Adjoint of control to initial state call gsi_bundleputvar ( wbundle_c, 't' , rv_tv, istatus ) call gsi_bundleputvar ( wbundle_c, 'ps', rv_ps, istatus ) + if (associated(rv_ext1)) & + call gsi_bundleputvar ( wbundle_c, 'ext1', rv_ext1, istatus ) + if (associated(rv_ext2)) & + call gsi_bundleputvar ( wbundle_c, 'ext2', rv_ext2, istatus ) ! call gsi_bundleputvar ( wbundle_c, 'q' , zero, istatus ) !$omp end parallel sections diff --git a/src/gsibec/gsi/get_gefs_ensperts_dualres.F90 b/src/gsibec/gsi/get_gefs_ensperts_dualres.F90 index cdd7590..141b294 100644 --- a/src/gsibec/gsi/get_gefs_ensperts_dualres.F90 +++ b/src/gsibec/gsi/get_gefs_ensperts_dualres.F90 @@ -99,7 +99,7 @@ subroutine get_gefs_ensperts_dualres (tau) ! integer(i_kind),dimension(grd_ens%nlat,grd_ens%nlon):: idum integer(i_kind) istatus,iret,i,ic2,ic3,j,k,n,mm1,iderivative,im,jm,km,m,ipic integer(i_kind) ipc3d(nc3d),ipc2d(nc2d) - integer(i_kind) ier + integer(i_kind) itv,iqv,ier ! integer(i_kind) il,jl logical ice,hydrometeor type(sub2grid_info) :: grd_tmp @@ -181,12 +181,18 @@ subroutine get_gefs_ensperts_dualres (tau) if (.not.q_hyb_ens) then !use RH call gsi_bundlegetpointer(en_read(n),'ps',ps,ier);istatus=ier - call gsi_bundlegetpointer(en_read(n),'t' ,tv,ier);istatus=istatus+ier - call gsi_bundlegetpointer(en_read(n),'q' ,q ,ier);istatus=istatus+ier + call gsi_bundlegetpointer(en_read(n),'t' ,tv,ier);itv=ier + call gsi_bundlegetpointer(en_read(n),'q' ,q ,ier);iqv=ier ! Compute RH ! Get 3d pressure field now on interfaces allocate(pri(im,jm,km+1)) - call general_getprs_glb(ps,tv,pri) + if(itv/=0) then + allocate(tv(im,jm,km)) + call general_getprs_glb(ps,tv,pri) + deallocate(tv) + else + call general_getprs_glb(ps,tv,pri) + endif allocate(prsl(im,jm,km),tsen(im,jm,km),qs(im,jm,km)) ! Get sensible temperature and 3d layer pressure if (idsl5 /= 2) then @@ -196,10 +202,19 @@ subroutine get_gefs_ensperts_dualres (tau) do i=1,im prsl(i,j,k)=((pri(i,j,k)**kap1-pri(i,j,k+1)**kap1)/& (kap1*(pri(i,j,k)-pri(i,j,k+1))))**kapr - tsen(i,j,k)= tv(i,j,k)/(one+fv*max(zero,q(i,j,k))) end do end do end do + if (itv==0.and.iqv==0) then +!$omp parallel do schedule(dynamic,1) private(k,j,i) + do k=1,km + do j=1,jm + do i=1,im + tsen(i,j,k)= tv(i,j,k)/(one+fv*max(zero,q(i,j,k))) + end do + end do + end do + end if else !$omp parallel do schedule(dynamic,1) private(k,j,i) do k=1,km diff --git a/src/gsibec/gsi/m_berror_stats.f90 b/src/gsibec/gsi/m_berror_stats.f90 index 6b0b9bf..838a233 100644 --- a/src/gsibec/gsi/m_berror_stats.f90 +++ b/src/gsibec/gsi/m_berror_stats.f90 @@ -585,7 +585,7 @@ subroutine bin_ call stop2(101) endif - write(6,*) myname_,'(PREWGT): read error amplitudes ', & + write(6,*) myname_,': read error amplitudes ', & '"',trim(berror_stats),'". ', & 'mype,nsigstat,nlatstat =', & mype,nsigstat,nlatstat @@ -838,15 +838,15 @@ subroutine setcoroz_(coroz,mype) if ( ierror/=0 ) return ! nothing to do ! sanity check - if ( mype==0 ) write(6,*) myname_,'(PREWGT): mype = ',mype + if ( mype==0 ) write(6,*) myname_,': mype = ',mype mlat=size(coroz,1) msig=size(coroz,2) if ( mlat/=nlat .or. msig/=nsig ) then - write(6,*) myname_,'(PREWGT): shape mismatching on PE ',mype - write(6,*) myname_,'(PREWGT): shape(coroz) = ',shape(coroz) - write(6,*) myname_,'(PREWGT): while expecting nlat = ',nlat - write(6,*) myname_,'(PREWGT): while expecting nsig = ',nsig + write(6,*) myname_,': shape mismatching on PE ',mype + write(6,*) myname_,': shape(coroz) = ',shape(coroz) + write(6,*) myname_,': while expecting nlat = ',nlat + write(6,*) myname_,': while expecting nsig = ',nsig call stop2(default_rc_) endif @@ -869,7 +869,7 @@ subroutine setcoroz_(coroz,mype) call mpi_allreduce(work_oz,work_oz1,(nsig+1)*npe,mpi_rtype,mpi_sum,& gsi_mpi_comm_world,ierror) if ( ierror/=0 ) then - write(6,*) myname_,'(PREWGT): MPI_allreduce() error on PE ',mype + write(6,*) myname_,': MPI_allreduce() error on PE ',mype call stop2(ierror) endif @@ -937,13 +937,13 @@ subroutine sethwlloz_(hwlloz,mype) real(r_kind) :: fact real(r_kind) :: s2u - if ( mype==0 ) write(6,*) myname_,'(PREWGT): mype = ',mype + if ( mype==0 ) write(6,*) myname_,': mype = ',mype s2u=(two*pi*rearth_equator)/nlon do k=1,nnnn1o k1=levs_id(k) if ( k1>0 ) then - if(mype==0) write(6,*) myname_,'(PREWGT): mype = ',mype, k1 + if(mype==0) write(6,*) myname_,': mype = ',mype, k1 if ( k1<=nsig*3/4 ) then ! fact=1./hwl fact=r40000/(r400*nlon) @@ -955,7 +955,7 @@ subroutine sethwlloz_(hwlloz,mype) endif enddo - if ( mype==0 ) write(6,*) myname_,'(PREWGT): mype = ',mype, 'finish sethwlloz_' + if ( mype==0 ) write(6,*) myname_,': mype = ',mype, 'finish sethwlloz_' return end subroutine sethwlloz_ @@ -1017,6 +1017,7 @@ subroutine setcorchem_(cname,corchem,rc) use guess_grids,only: ntguessig use guess_grids,only: ges_prsi ! interface pressures (kPa) + use gsi_metguess_mod, only: gsi_metguess_bundle use gsi_chemguess_mod,only: gsi_chemguess_bundle use gsi_bundlemod, only: gsi_bundlegetpointer @@ -1044,32 +1045,40 @@ subroutine setcorchem_(cname,corchem,rc) integer(i_kind) :: mlat,msig integer(i_kind) :: i,j,k,n,iptr,mm1 integer(i_kind) :: ierror + logical :: found_inchm, found_inmet rc=0 + found_inchm = .false. + found_inmet = .false. ! sanity check - if ( mype==0 ) write(6,*) myname_,'(PREWGT): mype = ',mype + if ( mype==0 ) write(6,*) myname_,': mype = ',mype ! Get information for how to use CO2 iptr=-1 if ( size(gsi_chemguess_bundle)>0 ) then ! check to see if bundle's allocated - call gsi_bundlegetpointer(gsi_chemguess_bundle(1),cname,iptr,ierror) + call gsi_bundlegetpointer(gsi_chemguess_bundle(1),cname,iptr,ierror) if ( ierror/=0 ) then rc=-2 ! field not found return endif + found_inchm=.true. else - rc=-1 ! chem not allocated - return + call gsi_bundlegetpointer(gsi_metguess_bundle(1),cname,iptr,ierror) + if ( ierror/=0 ) then + rc=-1 ! chem not allocated + return + endif + found_inmet=.true. endif mlat=size(corchem,1) msig=size(corchem,2) if ( mlat/=nlat .or. msig/=nsig ) then - write(6,*) myname_,'(PREWGT): shape mismatching on PE ',mype - write(6,*) myname_,'(PREWGT): shape(corchem',trim(cname),') = ',shape(corchem) - write(6,*) myname_,'(PREWGT): while expecting nlat = ',nlat - write(6,*) myname_,'(PREWGT): while expecting nsig = ',nsig + write(6,*) myname_,': shape mismatching on PE ',mype + write(6,*) myname_,': shape(corchem',trim(cname),') = ',shape(corchem) + write(6,*) myname_,': while expecting nlat = ',nlat + write(6,*) myname_,': while expecting nsig = ',nsig call stop2(default_rc_) endif @@ -1079,22 +1088,39 @@ subroutine setcorchem_(cname,corchem,rc) ! Calculate sums for constituent to estimate variance. mm1=mype+1 work_chem = zero - do k=1,nsig - do j=2,lon1+1 - do i=2,lat1+1 - work_chem(k,mm1) = work_chem(k,mm1) + gsi_chemguess_bundle(ntguessig)%r3(iptr)%q(i,j,k)* & - (ges_prsi(i,j,k,ntguessig)-ges_prsi(i,j,k+1,ntguessig)) - !_RT not sure yet how to handle scaling factor (rozcon) in general - !_RT rozcon*(ges_prsi(i,j,k,ntguessig)-ges_prsi(i,j,k+1,ntguessig)) + if ( found_inchm ) then + do k=1,nsig + do j=2,lon1+1 + do i=2,lat1+1 + work_chem(k,mm1) = work_chem(k,mm1) + gsi_chemguess_bundle(ntguessig)%r3(iptr)%q(i,j,k)* & + (ges_prsi(i,j,k,ntguessig)-ges_prsi(i,j,k+1,ntguessig)) + !_RT not sure yet how to handle scaling factor (rozcon) in general + !_RT rozcon*(ges_prsi(i,j,k,ntguessig)-ges_prsi(i,j,k+1,ntguessig)) + enddo enddo enddo - enddo + else if ( found_inmet ) then + do k=1,nsig + do j=2,lon1+1 + do i=2,lat1+1 + work_chem(k,mm1) = work_chem(k,mm1) + gsi_metguess_bundle(ntguessig)%r3(iptr)%q(i,j,k)* & + (ges_prsi(i,j,k,ntguessig)-ges_prsi(i,j,k+1,ntguessig)) + !_RT not sure yet how to handle scaling factor (rozcon) in general + !_RT rozcon*(ges_prsi(i,j,k,ntguessig)-ges_prsi(i,j,k+1,ntguessig)) + enddo + enddo + enddo + else + write(6,*) myname_,': should not be here ' + call stop2(ierror) + endif + work_chem(nsig+1,mm1)=float(lon1*lat1) call mpi_allreduce(work_chem,work_chem1,(nsig+1)*npe,mpi_rtype,mpi_sum,& gsi_mpi_comm_world,ierror) if ( ierror/=0 ) then - write(6,*) myname_,'(PREWGT): MPI_allreduce() error on PE ',mype + write(6,*) myname_,': MPI_allreduce() error on PE ',mype call stop2(ierror) endif @@ -1167,14 +1193,14 @@ subroutine sethwllchem_(hwll, mype) real(r_kind) :: s2u if (mype == 0) then - write(6,*) myname_, '(PREWGT): mype = ', mype + write(6,*) myname_, ': mype = ', mype end if s2u = (two*pi*rearth_equator)/nlon do k = 1,nnnn1o k1 = levs_id(k) if (k1 > 0) then - if (mype == 0) write(6,*) myname_, '(PREWGT): mype = ', mype, k1 + if (mype == 0) write(6,*) myname_, ': mype = ', mype, k1 ! make everything constant ! fact = real(k1,r_kind)**2._r_kind fact = 1._r_kind @@ -1184,7 +1210,7 @@ subroutine sethwllchem_(hwll, mype) end do if (mype == 0) then - write(6,*) myname_, '(PREWGT): mype = ', mype, 'finish sethwllchem_' + write(6,*) myname_, ': mype = ', mype, 'finish sethwllchem_' end if end subroutine sethwllchem_ diff --git a/src/gsibec/gsigeos/m_nc_GEOSens.f90 b/src/gsibec/gsigeos/m_nc_GEOSens.f90 index 760726e..946e9ad 100644 --- a/src/gsibec/gsigeos/m_nc_GEOSens.f90 +++ b/src/gsibec/gsigeos/m_nc_GEOSens.f90 @@ -1,8 +1,10 @@ module m_nc_GEOSens use netcdf +use mpeu_util, only: getindex implicit none private +public :: nc_GEOSens_vars_set public :: nc_GEOSens_vars_init public :: nc_GEOSens_vars_final public :: nc_GEOSens_vars_comp @@ -16,18 +18,24 @@ module m_nc_GEOSens public :: nc_GEOSens_geos2gsi public :: nc_GEOSens_gsi2geos +! The following type is essentially the CV but the varnames +! it operates on a those in file, as opposed to the GSI +! names - see anavinfo. Because this is controlled in +! anavinfo, the ordering of the fields is that in the CV. +! See also, set_ routine below and its input fields names +! arrays. type nc_GEOSens_vars logical :: initialized=.false. integer :: nlon,nlat,nsig logical :: gsiset=.false. real(4),pointer,dimension(:):: ak,bk - real(4),pointer,dimension(:,:,:):: dp - real(4),pointer,dimension(:,:,:):: tv - real(4),pointer,dimension(:,:,:):: u,v - real(4),pointer,dimension(:,:,:):: qv - real(4),pointer,dimension(:,:,:):: qi,ql,qr,qs - real(4),pointer,dimension(:,:,:):: oz - real(4),pointer,dimension(:,:) :: ps,ts + real(4),pointer,dimension(:,:,:,:):: ptr3d + real(4),pointer,dimension(:,:,:) :: ptr2d +! + integer :: nv2d = -1 + integer :: nv3d = -1 + character(len=5),allocatable :: fvars2d(:) + character(len=5),allocatable :: fvars3d(:) ! real(4),pointer,dimension(:) :: v1d real(4),pointer,dimension(:,:) :: v2d @@ -39,16 +47,6 @@ module m_nc_GEOSens real, parameter:: mbar_per_Pa = 0.01 ! mb to Pa real, parameter:: Pa_per_kPa = 1000.0 -integer, parameter :: nv2d = 2 -character(len=4),parameter :: cvars2d(nv2d) = (/ 'ps ', 'ts ' /) - -integer, parameter :: nv3d = 10 -character(len=5),parameter :: cvars3d(nv3d) = (/ & - 'tv ', 'u ', 'v ', & - 'sphu ', 'qitot', 'qltot', & - 'qrtot', 'qstot', 'ozone', & - 'delp '& - /) interface nc_GEOSens_dims; module procedure & read_dims_ ; end interface @@ -56,6 +54,8 @@ module m_nc_GEOSens read_GEOSens_ ; end interface interface nc_GEOSens_write; module procedure & write_GEOSens_ ; end interface +interface nc_GEOSens_vars_set; module procedure & + set_vars_ ; end interface interface nc_GEOSens_vars_init; module procedure & init_GEOSens_vars_ ; end interface interface nc_GEOSens_vars_final; module procedure & @@ -82,6 +82,24 @@ module m_nc_GEOSens end interface contains +subroutine set_vars_(fvars2d,fvars3d,bvars) +implicit none +character(*), intent(in) :: fvars2d(:) ! in-file-varname for 2d-fiels per anavinfo +character(*), intent(in) :: fvars3d(:) ! in-file-varname for 3d-fiels per anavinfo +type(nc_GEOSens_vars),intent(inout) :: bvars + +bvars%nv2d = size(fvars2d) +bvars%nv3d = size(fvars3d) +if (bvars%nv2d>0) then + if(.not.allocated(bvars%fvars2d)) allocate(bvars%fvars2d(bvars%nv2d)) + bvars%fvars2d = fvars2d +endif +if (bvars%nv3d>0) then + if(.not.allocated(bvars%fvars3d)) allocate(bvars%fvars3d(bvars%nv2d)) + bvars%fvars3d = fvars3d +endif +end subroutine set_vars_ + subroutine read_dims_ (fname,nlat,nlon,nlev,rc, myid,root) implicit none character(len=*), intent(in) :: fname ! input filename @@ -209,90 +227,28 @@ subroutine read_GEOSens_ (fname,bvars,rc, myid,root, gsiset) ! Read data to file allocate(data_in(nlon,nlat,1)) - do nv = 1, nv2d - call check_( nf90_inq_varid(ncid, trim(cvars2d(nv)), varid), rc, mype_, root_ ) + do nv = 1, bvars%nv2d + call check_( nf90_inq_varid(ncid, trim(bvars%fvars2d(nv)), varid), rc, mype_, root_ ) call check_( nf90_get_var(ncid, varid, data_in(:,:,1)), rc, mype_, root_ ) if (bvars%gsiset) then - if(trim(cvars2d(nv))=="ps" ) bvars%ps = transpose(data_in(:,:,1)) - if(trim(cvars2d(nv))=="ts" ) bvars%ts = transpose(data_in(:,:,1)) + bvars%ptr2d(:,:,nv) = transpose(data_in(:,:,1)) else - if(trim(cvars2d(nv))=="ps" ) bvars%ps = data_in(:,:,1) - if(trim(cvars2d(nv))=="ts" ) bvars%ts = data_in(:,:,1) + bvars%ptr2d(:,:,nv) = data_in(:,:,1) endif enddo deallocate(data_in) ! allocate(data_in(nlon,nlat,nlev)) - do nv = 1, nv3d - call check_( nf90_inq_varid(ncid, trim(cvars3d(nv)), varid), rc, mype_, root_ ) + do nv = 1, bvars%nv3d + call check_( nf90_inq_varid(ncid, trim(bvars%fvars3d(nv)), varid), rc, mype_, root_ ) call check_( nf90_get_var(ncid, varid, data_in(:,:,:)), rc, mype_, root_ ) if (bvars%gsiset) then - if(trim(cvars3d(nv))=="delp") then - do kk=1,bvars%nsig - bvars%dp(:,:,kk) = transpose(data_in(:,:,kk)) - enddo - endif - if(trim(cvars3d(nv))=="tv" ) then - do kk=1,bvars%nsig - bvars%tv(:,:,kk) = transpose(data_in(:,:,kk)) - enddo - endif - if(trim(cvars3d(nv))=="u" ) then - do kk=1,bvars%nsig - bvars%u(:,:,kk) = transpose(data_in(:,:,kk)) - enddo - endif - if(trim(cvars3d(nv))=="v" ) then - do kk=1,bvars%nsig - bvars%v(:,:,kk) = transpose(data_in(:,:,kk)) - enddo - endif -! - if(trim(cvars3d(nv))=="sphu" ) then - do kk=1,bvars%nsig - bvars%qv(:,:,kk) = transpose(data_in(:,:,kk)) - enddo - endif - if(trim(cvars3d(nv))=="qitot") then - do kk=1,bvars%nsig - bvars%qi(:,:,kk) = transpose(data_in(:,:,kk)) - enddo - endif - if(trim(cvars3d(nv))=="qltot") then - do kk=1,bvars%nsig - bvars%ql(:,:,kk) = transpose(data_in(:,:,kk)) - enddo - endif - if(trim(cvars3d(nv))=="qrtot") then - do kk=1,bvars%nsig - bvars%qr(:,:,kk) = transpose(data_in(:,:,kk)) - enddo - endif - if(trim(cvars3d(nv))=="qstot") then - do kk=1,bvars%nsig - bvars%qs(:,:,kk) = transpose(data_in(:,:,kk)) - enddo - endif -! - if(trim(cvars3d(nv))=="ozone") then - do kk=1,bvars%nsig - bvars%oz(:,:,kk) = transpose(data_in(:,:,kk)) - enddo - endif + do kk=1,nlev + bvars%ptr3d(:,:,kk,nv) = transpose(data_in(:,:,kk)) + enddo else - if(trim(cvars3d(nv))=="delp") bvars%dp = data_in(:,:,:) - if(trim(cvars3d(nv))=="tv" ) bvars%tv = data_in(:,:,:) - if(trim(cvars3d(nv))=="u" ) bvars%u = data_in(:,:,:) - if(trim(cvars3d(nv))=="v" ) bvars%v = data_in(:,:,:) -! - if(trim(cvars3d(nv))=="sphu" ) bvars%qv = data_in(:,:,:) - if(trim(cvars3d(nv))=="qitot") bvars%qi = data_in(:,:,:) - if(trim(cvars3d(nv))=="qltot") bvars%ql = data_in(:,:,:) - if(trim(cvars3d(nv))=="qrtot") bvars%qr = data_in(:,:,:) - if(trim(cvars3d(nv))=="qstot") bvars%qs = data_in(:,:,:) -! - if(trim(cvars3d(nv))=="ozone") bvars%oz = data_in(:,:,:) + bvars%ptr3d(:,:,:,nv) = data_in(:,:,:) endif ! enddo @@ -392,13 +348,13 @@ subroutine write_GEOSens_ (fname,bvars,lats,lons,rc, myid,root,plevs) dimids = (/ x_dimid, y_dimid, z_dimid /) ! Define variables. - allocate(varid2d(nv2d)) - do nv = 1, nv2d - call check_( nf90_def_var(ncid, trim(cvars2d(nv)), NF90_REAL, (/ x_dimid, y_dimid /), varid2d(nv)), rc, mype_, root_ ) + allocate(varid2d(bvars%nv2d)) + do nv = 1, bvars%nv2d + call check_( nf90_def_var(ncid, trim(bvars%fvars2d(nv)), NF90_REAL, (/ x_dimid, y_dimid /), varid2d(nv)), rc, mype_, root_ ) enddo - allocate(varid3d(nv3d)) - do nv = 1, nv3d - call check_( nf90_def_var(ncid, trim(cvars3d(nv)), NF90_REAL, (/ x_dimid, y_dimid, z_dimid /), varid3d(nv)), rc, mype_, root_ ) + allocate(varid3d(bvars%nv3d)) + do nv = 1, bvars%nv3d + call check_( nf90_def_var(ncid, trim(bvars%fvars3d(nv)), NF90_REAL, (/ x_dimid, y_dimid, z_dimid /), varid3d(nv)), rc, mype_, root_ ) enddo ! End define mode. This tells netCDF we are done defining metadata. @@ -420,26 +376,14 @@ subroutine write_GEOSens_ (fname,bvars,lats,lons,rc, myid,root,plevs) ! Write data to file allocate(data_out(nlon,nlat,1)) - do nv = 1, nv2d - if(trim(cvars2d(nv))=="ps" ) data_out(:,:,1) = bvars%ps - if(trim(cvars2d(nv))=="ts" ) data_out(:,:,1) = bvars%ts + do nv = 1, bvars%nv2d + data_out(:,:,1) = bvars%ptr2d(:,:,nv) call check_( nf90_put_var(ncid, varid2d(nv), data_out(:,:,1)), rc, mype_, root_) enddo deallocate(data_out) allocate(data_out(nlon,nlat,nlev)) - do nv = 1, nv3d - if(trim(cvars3d(nv))=="delp") data_out(:,:,:) = bvars%dp - if(trim(cvars3d(nv))=="tv" ) data_out(:,:,:) = bvars%tv - if(trim(cvars3d(nv))=="u" ) data_out(:,:,:) = bvars%u - if(trim(cvars3d(nv))=="v" ) data_out(:,:,:) = bvars%v -! - if(trim(cvars2d(nv))=="sphu" ) data_out(:,:,:) = bvars%qv - if(trim(cvars2d(nv))=="qitot") data_out(:,:,:) = bvars%qi - if(trim(cvars2d(nv))=="qltot") data_out(:,:,:) = bvars%ql - if(trim(cvars2d(nv))=="qrtot") data_out(:,:,:) = bvars%qr - if(trim(cvars2d(nv))=="qstot") data_out(:,:,:) = bvars%qs -! - if(trim(cvars2d(nv))=="ozone") data_out(:,:,:) = bvars%oz + do nv = 1, bvars%nv3d + data_out(:,:,:) = bvars%ptr3d(:,:,:,nv) ! call check_( nf90_put_var(ncid, varid3d(nv), data_out(:,:,:)), rc, mype_, root_ ) enddo @@ -463,6 +407,8 @@ subroutine init_GEOSens_vars_(vr,nlon,nlat,nsig,gsi) type(nc_GEOSens_vars) vr logical,intent(in),optional :: gsi + integer :: ii + if(vr%initialized) return if(present(gsi)) then @@ -474,16 +420,13 @@ subroutine init_GEOSens_vars_(vr,nlon,nlat,nsig,gsi) vr%nsig=nsig ! allocate single precision arrays +! why the if? if (vr%gsiset) then - allocate(vr%tv(nlat,nlon,nsig),vr%u (nlat,nlon,nsig),vr%v (nlat,nlon,nsig),vr%qv(nlat,nlon,nsig),& - vr%qi(nlat,nlon,nsig),vr%ql(nlat,nlon,nsig),vr%qr(nlat,nlon,nsig),vr%qs(nlat,nlon,nsig),& - vr%oz(nlat,nlon,nsig),vr%dp(nlat,nlon,nsig) ) - allocate(vr%ps(nlat,nlon),vr%ts(nlat,nlon)) + if(vr%nv3d>0) allocate(vr%ptr3d(nlat,nlon,nsig,vr%nv3d)) + if(vr%nv2d>0) allocate(vr%ptr2d(nlat,nlon,vr%nv2d)) else - allocate(vr%tv(nlon,nlat,nsig),vr%u (nlon,nlat,nsig),vr%v (nlon,nlat,nsig),vr%qv(nlon,nlat,nsig),& - vr%qi(nlon,nlat,nsig),vr%ql(nlon,nlat,nsig),vr%qr(nlon,nlat,nsig),vr%qs(nlon,nlat,nsig),& - vr%oz(nlon,nlat,nsig),vr%dp(nlon,nlat,nsig) ) - allocate(vr%ps(nlon,nlat),vr%ts(nlon,nlat)) + if(vr%nv3d>0) allocate(vr%ptr3d(nlat,nlon,nsig,vr%nv3d)) + if(vr%nv2d>0) allocate(vr%ptr2d(nlat,nlon,vr%nv2d)) endif vr%initialized=.true. end subroutine init_GEOSens_vars_ @@ -492,10 +435,8 @@ subroutine final_GEOSens_vars_(vr) type(nc_GEOSens_vars) vr ! deallocate arrays if(.not. vr%initialized) return - deallocate(vr%tv,vr%u ,vr%v ,vr%qv, & - vr%qi,vr%ql,vr%qr,vr%qs,& - vr%oz,vr%dp) - deallocate(vr%ps,vr%ts) + if(vr%nv2d>0) deallocate(vr%ptr2d) + if(vr%nv3d>0) deallocate(vr%ptr3d) vr%initialized=.false. end subroutine final_GEOSens_vars_ @@ -505,7 +446,7 @@ subroutine comp_GEOSens_vars_(va,vb,rc, myid,root) integer, intent(out) :: rc integer, intent(in), optional :: myid,root ! accommodate MPI calling programs character(len=*), parameter :: myname_ = myname//"::comp_GEOSens_vars_" - integer :: ii,jj + integer :: ii,jj,nv logical :: verbose, failed real :: tolerance = 10.e-10 integer, allocatable :: ier(:) @@ -527,20 +468,26 @@ subroutine comp_GEOSens_vars_(va,vb,rc, myid,root) return endif - allocate(ier(nv2d+nv3d)) - ii=0;ier=0 - ii=ii+1; if(abs(sum(va%dp - vb%dp)) >tolerance) ier(ii)=ii - ii=ii+1; if(abs(sum(va%tv - vb%tv)) >tolerance) ier(ii)=ii - ii=ii+1; if(abs(sum(va%u - vb%u )) >tolerance) ier(ii)=ii - ii=ii+1; if(abs(sum(va%v - vb%v )) >tolerance) ier(ii)=ii - ii=ii+1; if(abs(sum(va%qv - vb%qv)) >tolerance) ier(ii)=ii - ii=ii+1; if(abs(sum(va%qi - vb%qi)) >tolerance) ier(ii)=ii - ii=ii+1; if(abs(sum(va%ql - vb%ql)) >tolerance) ier(ii)=ii - ii=ii+1; if(abs(sum(va%qr - vb%qr)) >tolerance) ier(ii)=ii - ii=ii+1; if(abs(sum(va%qs - vb%qs)) >tolerance) ier(ii)=ii - ii=ii+1; if(abs(sum(va%oz - vb%oz)) >tolerance) ier(ii)=ii - ii=ii+1; if(abs(sum(va%ps - vb%ps)) >tolerance) ier(ii)=ii - ii=ii+1; if(abs(sum(va%ts - vb%ts)) >tolerance) ier(ii)=ii + if(va%nv2d /= vb%nv2d .or. va%nv3d /= vb%nv3d) then + rc=1 + if(myid==root) then + print *, 'nv2d(va) = ', va%nv2d, 'nv2d(vb) = ', vb%nv2d + print *, 'nv3d(va) = ', va%nv3d, 'nv3d(vb) = ', vb%nv3d + print *, myname_, 'Inconsistent dimensions, aborting ... ' + endif + endif + allocate(ier(va%nv2d+va%nv3d)) + ier=0 + + ii=0 + do nv=1,va%nv3d + ii=ii+1 + if(abs(sum(va%ptr3d(:,:,:,nv) - vb%ptr3d(:,:,:,nv))) >tolerance) ier(ii)=ii + enddo + do nv=1,va%nv2d + ii=ii+1 + if(abs(sum(va%ptr2d(:,:,nv) - vb%ptr2d(:,:,nv))) >tolerance) ier(ii)=ii + enddo failed=.false. do jj=1,ii if(ier(jj)/=0.and.verbose) then @@ -558,7 +505,7 @@ subroutine copy_(ivars,ovars,rc) type(nc_GEOSens_vars) ivars type(nc_GEOSens_vars) ovars integer, intent(out) :: rc - integer :: kk + integer :: kk,nv rc=0 if (ovars%nlon/=ivars%nlon .or. & @@ -570,34 +517,21 @@ subroutine copy_(ivars,ovars,rc) endif if (ivars%gsiset .neqv. ovars%gsiset ) then - do kk=1,ovars%nsig - ovars%dp(:,:,kk) = transpose(ivars%dp(:,:,kk)) - ovars%tv(:,:,kk) = transpose(ivars%tv(:,:,kk)) - ovars%u (:,:,kk) = transpose(ivars%u (:,:,kk)) - ovars%v (:,:,kk) = transpose(ivars%v (:,:,kk)) - ovars%qv(:,:,kk) = transpose(ivars%qv(:,:,kk)) - ovars%qi(:,:,kk) = transpose(ivars%qi(:,:,kk)) - ovars%ql(:,:,kk) = transpose(ivars%ql(:,:,kk)) - ovars%qr(:,:,kk) = transpose(ivars%qr(:,:,kk)) - ovars%qs(:,:,kk) = transpose(ivars%qs(:,:,kk)) - ovars%oz(:,:,kk) = transpose(ivars%oz(:,:,kk)) + do nv=1,ovars%nv3d + do kk=1,ovars%nsig + ovars%ptr3d(:,:,kk,nv) = transpose(ivars%ptr3d(:,:,kk,nv)) + enddo + enddo + do nv=1,ovars%nv2d + ovars%ptr2d(:,:,nv) = transpose(ivars%ptr2d(:,:,nv)) enddo - ovars%ps = transpose(ivars%ps) - ovars%ts = transpose(ivars%ts) else - ovars%dp = ivars%dp - ovars%tv = ivars%tv - ovars%u = ivars%u - ovars%v = ivars%v - ovars%qv = ivars%qv - ovars%qi = ivars%qi - ovars%ql = ivars%ql - ovars%qr = ivars%qr - ovars%qs = ivars%qs - ovars%oz = ivars%oz - - ovars%ps = ivars%ps - ovars%ts = ivars%ts + do nv=1,ovars%nv3d + ovars%ptr3d(:,:,:,nv) = ivars%ptr3d(:,:,:,nv) + enddo + do nv=1,ovars%nv2d + ovars%ptr2d(:,:,nv) = ivars%ptr2d(:,:,nv) + enddo endif end subroutine copy_ @@ -608,13 +542,11 @@ subroutine get_pointer_2d_ (vname, bvars, ptr, rc ) type(nc_GEOSens_vars) bvars real(4),pointer,intent(inout) :: ptr(:,:) integer,intent(out) :: rc +integer :: id rc=-1 -if(trim(vname)=='ps') then - ptr => bvars%ps - rc=0 -endif -if(trim(vname)=='ts') then - ptr => bvars%ts +id=getindex(bvars%fvars2d,trim(vname)) +if (id>0) then + ptr = bvars%ptr2d(:,:,id) rc=0 endif end subroutine get_pointer_2d_ @@ -625,77 +557,12 @@ subroutine get_pointer_3d_ (vname, bvars, ptr, rc ) type(nc_GEOSens_vars) bvars real(4),pointer,intent(inout) :: ptr(:,:,:) integer,intent(out) :: rc -character(len=5) :: var +integer :: id rc=-1 -! -var='delp' -if(trim(vname)==trim(var)) then - ptr => bvars%dp - rc=0 - return -endif -! -var='tv' -if(trim(vname)==trim(var)) then - ptr => bvars%tv - rc=0 - return -endif -! -var='u' -if(trim(vname)==trim(var)) then - ptr => bvars%u - rc=0 - return -endif -! -var='v' -if(trim(vname)==trim(var)) then - ptr => bvars%v - rc=0 - return -endif -! -var='sphu' -if(trim(vname)==trim(var)) then - ptr => bvars%qv - rc=0 - return -endif -! -var='qitot' -if(trim(vname)==trim(var)) then - ptr => bvars%qi - rc=0 - return -endif -! -var='qltot' -if(trim(vname)==trim(var)) then - ptr => bvars%ql - rc=0 - return -endif -! -var='qrtot' -if(trim(vname)==trim(var)) then - ptr => bvars%qr +id=getindex(bvars%fvars3d,trim(vname)) +if (id>0) then + ptr = bvars%ptr3d(:,:,:,id) rc=0 - return -endif -! -var='qstot' -if(trim(vname)==trim(var)) then - ptr => bvars%qs - rc=0 - return -endif -! -var='ozone' -if(trim(vname)==trim(var)) then - ptr => bvars%oz - rc=0 - return endif end subroutine get_pointer_3d_ @@ -713,12 +580,28 @@ end subroutine check_ subroutine geos2gsi_ (x) implicit none type(nc_GEOSens_vars) x + integer id !==> input ps in Pa - convert to hPa(mb) ! x%grid%ak = x%grid%ak / Pa_per_kPa - x%ps = x%ps / Pa_per_kPa - x%dp = x%dp / Pa_per_kPa - x%oz = x%oz / PPMV2GpG + id=getindex(x%fvars2d,'ps') + if (id>0) then + x%ptr2d(:,:,id) = x%ptr2d(:,:,id) / Pa_per_kPa + endif + id=getindex(x%fvars3d,'delp') + if (id>0) then + x%ptr3d(:,:,:,id) = x%ptr3d(:,:,:,id) / Pa_per_kPa + else + id=getindex(x%fvars3d,'dp') + if(id>0) x%ptr3d(:,:,:,id) = x%ptr3d(:,:,:,id) / Pa_per_kPa + endif + id=getindex(x%fvars3d,'ozone') + if (id>0) then + x%ptr3d(:,:,:,id) = x%ptr3d(:,:,:,id) / PPMV2GpG + else + id=getindex(x%fvars3d,'oz') + if(id>0) x%ptr3d(:,:,:,id) = x%ptr3d(:,:,:,id) / PPMV2GpG + endif ! need flip so localization function applies equaly to EnKF and Hybrid-GSI call flip_(x) ! still need a transpose @@ -727,12 +610,28 @@ end subroutine geos2gsi_ subroutine gsi2geos_ (x) implicit none type(nc_GEOSens_vars) x + integer :: id !==> input ps in mbar - convert to Pa ! x%grid%ak = x%grid%ak * Pa_per_kPa - x%ps = x%ps * Pa_per_kPa - x%dp = x%dp * Pa_per_kPa - x%oz = x%oz * PPMV2GpG + id=getindex(x%fvars2d,'ps') + if (id>0) then + x%ptr2d(:,:,id) = x%ptr2d(:,:,id) * Pa_per_kPa + endif + id=getindex(x%fvars3d,'delp') + if (id>0) then + x%ptr3d(:,:,:,id) = x%ptr3d(:,:,:,id) * Pa_per_kPa + else + id=getindex(x%fvars3d,'dp') + if (id>0) x%ptr3d(:,:,:,id) = x%ptr3d(:,:,:,id) * Pa_per_kPa + endif + id=getindex(x%fvars3d,'ozone') + if (id>0) then + x%ptr3d(:,:,:,id) = x%ptr3d(:,:,:,id) * PPMV2GpG + else + id=getindex(x%fvars3d,'oz') + if(id>0) x%ptr3d(:,:,:,id) = x%ptr3d(:,:,:,id) * PPMV2GpG + endif ! need flip so localization function applies equaly to EnKF and Hybrid-GSI call flip_(x) ! still need a transpose @@ -741,43 +640,19 @@ end subroutine gsi2geos_ subroutine flip_(x) implicit none type(nc_GEOSens_vars) x - integer im,jm,km + integer im,jm,km,nv im=x%nlon jm=x%nlat km=x%nsig ! - call hflip2_(x%ps,im,jm,x%gsiset) - call hflip2_(x%ts,im,jm,x%gsiset) + do nv=1,x%nv2d + call hflip2_(x%ptr2d(:,:,nv),im,jm,x%gsiset) + enddo ! - call hflip3_(x%dp,im,jm,km,x%gsiset) - call vflip_ (x%dp,im,jm,km) - - call hflip3_(x%tv,im,jm,km,x%gsiset) - call vflip_ (x%tv,im,jm,km) - - call hflip3_(x%u ,im,jm,km,x%gsiset) - call vflip_ (x%u ,im,jm,km) - - call hflip3_(x%v ,im,jm,km,x%gsiset) - call vflip_ (x%v ,im,jm,km) - - call hflip3_(x%qv,im,jm,km,x%gsiset) - call vflip_ (x%qv,im,jm,km) - - call hflip3_(x%qi,im,jm,km,x%gsiset) - call vflip_ (x%qi,im,jm,km) - - call hflip3_(x%ql,im,jm,km,x%gsiset) - call vflip_ (x%ql,im,jm,km) - - call hflip3_(x%qr,im,jm,km,x%gsiset) - call vflip_ (x%qr,im,jm,km) - - call hflip3_(x%qs,im,jm,km,x%gsiset) - call vflip_ (x%qs,im,jm,km) - - call hflip3_(x%oz,im,jm,km,x%gsiset) - call vflip_ (x%oz,im,jm,km) + do nv=1,x%nv3d + call hflip3_(x%ptr3d(:,:,:,nv),im,jm,km,x%gsiset) + call vflip_ (x%ptr3d(:,:,:,nv),im,jm,km) + enddo end subroutine flip_ subroutine hflip3_ ( q,im,jm,km, gsi ) @@ -858,25 +733,27 @@ subroutine summary_(x,myid) implicit none type(nc_GEOSens_vars) x integer, intent(in), optional :: myid - integer lu,nxy,nxyz + integer lu,nxy,nxyz,nv nxy = x%nlat*x%nlon nxyz = nxy*x%nsig lu=6 if(present(myid)) lu=myid - write(lu,'(a)') "================================================" - write(lu,'(a)') "var min max mean stddev " - write(lu,'(a)') "================================================" - write(lu,'(a4,1p,4(e10.3,1x))') 'ps', minval(x%ps), maxval(x%ps), sum(x%ps)/nxy, stddev_(x%ps) - write(lu,'(a4,1p,4(e10.3,1x))') 'ts', minval(x%ts), maxval(x%ts), sum(x%ts)/nxy, stddev_(x%ts) - write(lu,'(a4,1p,4(e10.3,1x))') 'dp', minval(x%dp), maxval(x%dp), sum(x%dp)/nxyz, stddev_(x%dp) - write(lu,'(a4,1p,4(e10.3,1x))') 'tv', minval(x%tv), maxval(x%tv), sum(x%tv)/nxyz, stddev_(x%tv) - write(lu,'(a4,1p,4(e10.3,1x))') 'qv', minval(x%qv), maxval(x%qv), sum(x%qv)/nxyz, stddev_(x%qv) - write(lu,'(a4,1p,4(e10.3,1x))') 'qi', minval(x%qi), maxval(x%qi), sum(x%qi)/nxyz, stddev_(x%qi) - write(lu,'(a4,1p,4(e10.3,1x))') 'ql', minval(x%ql), maxval(x%ql), sum(x%ql)/nxyz, stddev_(x%ql) - write(lu,'(a4,1p,4(e10.3,1x))') 'qr', minval(x%qr), maxval(x%qr), sum(x%qr)/nxyz, stddev_(x%qr) - write(lu,'(a4,1p,4(e10.3,1x))') 'qs', minval(x%qs), maxval(x%qs), sum(x%qs)/nxyz, stddev_(x%qs) - write(lu,'(a4,1p,4(e10.3,1x))') 'oz', minval(x%oz), maxval(x%oz), sum(x%oz)/nxyz, stddev_(x%oz) - write(lu,'(a)') "================================================" + write(lu,'(a)') "=================================================" + write(lu,'(a)') "var min max mean stddev " + write(lu,'(a)') "=================================================" + do nv=1,x%nv2d + write(lu,'(a5,1p,4(e10.3,1x))') trim(x%fvars2d(nv)), minval(x%ptr2d(:,:,nv)), & + maxval(x%ptr2d(:,:,nv)), & + sum(x%ptr2d(:,:,nv))/nxy, & + stddev_(x%ptr2d(:,:,nv)) + enddo + do nv=1,x%nv3d + write(lu,'(a5,1p,4(e10.3,1x))') trim(x%fvars3d(nv)), minval(x%ptr3d(:,:,:,nv)), & + maxval(x%ptr3d(:,:,:,nv)), & + sum(x%ptr3d(:,:,:,nv))/nxy, & + stddev_(x%ptr3d(:,:,:,nv)) + enddo + write(lu,'(a)') "=================================================" end subroutine summary_ real function stddev2_(x) diff --git a/src/gsibec/gsigeos/m_read_geosens.F90 b/src/gsibec/gsigeos/m_read_geosens.F90 index 18d4caf..5e53544 100644 --- a/src/gsibec/gsigeos/m_read_geosens.F90 +++ b/src/gsibec/gsigeos/m_read_geosens.F90 @@ -13,6 +13,7 @@ module m_read_geosens use gridmod,only: nsig use control_vectors, only: nc2d,nc3d use control_vectors, only: cvars2d,cvars3d + use control_vectors, only: fvars2d,fvars3d use general_sub2grid_mod, only: sub2grid_info use m_grid2sub1var, only: grid2sub1var @@ -34,6 +35,7 @@ module m_read_geosens contains subroutine read_geosens_(sgrid,xx,filename,npe,mype,root,nreaders) + use m_nc_GEOSens, only: nc_GEOSens_vars_set use m_nc_GEOSens, only: nc_GEOSens_vars use m_nc_GEOSens, only: nc_GEOSens_read use m_nc_GEOSens, only: nc_GEOSens_vars_final @@ -67,6 +69,7 @@ subroutine read_geosens_(sgrid,xx,filename,npe,mype,root,nreaders) nlon=sgrid%nlon iglobal=sgrid%iglobal + call nc_GEOSens_vars_set(fvars2d,fvars3d,evars) ! wired for test call init_() mm1=mype+1 @@ -109,6 +112,7 @@ subroutine read_geosens_(sgrid,xx,filename,npe,mype,root,nreaders) call final_() + return contains @@ -166,137 +170,27 @@ subroutine read_(n,proc1) ke=0 do nv=1,nc3d - if (cvars3d(nv)=='sf') then - do k=1,nsig - ii =0 - do j=1,nlon - do i=1,nlat - ii = ii +1 - z4all(ii,ke+k)=evars%u(i,j,k) - enddo - enddo - enddo - ke=ke+nsig - endif - if (cvars3d(nv)=='vp') then - do k=1,nsig - ii =0 - do j=1,nlon - do i=1,nlat - ii = ii +1 - z4all(ii,ke+k)=evars%v(i,j,k) - enddo - enddo - enddo - ke=ke+nsig - endif - if (cvars3d(nv)=='t') then - do k=1,nsig - ii =0 - do j=1,nlon - do i=1,nlat - ii = ii +1 - z4all(ii,ke+k)=evars%tv(i,j,k) - enddo - enddo - enddo - ke=ke+nsig - endif - if (cvars3d(nv)=='q') then - do k=1,nsig - ii =0 - do j=1,nlon - do i=1,nlat - ii = ii +1 - z4all(ii,ke+k)=evars%qv(i,j,k) - enddo - enddo - enddo - ke=ke+nsig - endif - if (cvars3d(nv)=='oz') then - do k=1,nsig - ii =0 - do j=1,nlon - do i=1,nlat - ii = ii +1 - z4all(ii,ke+k)=evars%oz(i,j,k) - enddo - enddo - enddo - ke=ke+nsig - endif - if (cvars3d(nv)=='qi') then - do k=1,nsig - ii =0 - do j=1,nlon - do i=1,nlat - ii = ii +1 - z4all(ii,ke+k)=evars%qi(i,j,k) - enddo - enddo - enddo - ke=ke+nsig - endif - if (cvars3d(nv)=='ql') then - do k=1,nsig - ii =0 - do j=1,nlon - do i=1,nlat - ii = ii +1 - z4all(ii,ke+k)=evars%ql(i,j,k) - enddo - enddo - enddo - ke=ke+nsig - endif - if (cvars3d(nv)=='qr') then - do k=1,nsig - ii =0 - do j=1,nlon - do i=1,nlat - ii = ii +1 - z4all(ii,ke+k)=evars%qr(i,j,k) - enddo - enddo + do k=1,nsig + ii =0 + do j=1,nlon + do i=1,nlat + ii = ii +1 + z4all(ii,ke+k)=evars%ptr3d(i,j,k,nv) enddo - ke=ke+nsig - endif - if (cvars3d(nv)=='qs') then - do k=1,nsig - ii =0 - do j=1,nlon - do i=1,nlat - ii = ii +1 - z4all(ii,ke+k)=evars%qs(i,j,k) - end do - end do enddo - ke=ke+nsig - endif + enddo + ke=ke+nsig enddo do nv=1,nc2d - if (cvars2d(nv)=='ps') then - ii=0 - do j=1,nlon - do i=1,nlat - ii = ii+1 - z4all(ii,ke+1)=evars%ps(i,j) - enddo - enddo - ke=ke+1 - endif - if (cvars2d(nv)=='ts') then - ii=0 - do j=1,nlon - do i=1,nlat - ii = ii+1 - z4all(ii,ke+1)=evars%ts(i,j) - enddo - enddo - ke=ke+1 - endif + ii=0 + do j=1,nlon + do i=1,nlat + ii = ii+1 + z4all(ii,ke+1)=evars%ptr2d(i,j,nv) + enddo + enddo + ke=ke+1 enddo call nc_GEOSens_vars_final(evars)