Feature/rtodling/add extinctions

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

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
 

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

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 )

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
 

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