Merge pull request NOAA-EMC#4 from guoqing-noaa/merge

update setupt.f90, update_guess.f90 etc based on latest GSL 2mDA code

src/gsi/gsd_update_mod.f90

@@ -15,7 +15,7 @@ module gsd_update_mod
 ! subroutines included:
 !   sub gsd_update_soil_tq  - change surface and soil based on analysis increment
 !   sub gsd_limit_ocean_q   - limits to analysis increments over oceans 
-!   sub gsd_update_th2      - adjust 2-m t based on analysis increment
+!   sub gsd_update_t2m      - adjust 2-m t based on analysis increment
 !   sub gsd_update_q2       - adjust 2-m q based on analysis increment
 !
 ! Variable Definitions:
@@ -30,7 +30,7 @@ module gsd_update_mod
 ! set subroutines to public
   public :: gsd_update_soil_tq
   public :: gsd_limit_ocean_q
-  public :: gsd_update_th2
+  public :: gsd_update_t2m
   public :: gsd_update_q2
   public :: gsd_gen_coast_prox
 ! set passed variables to public
@@ -487,10 +487,10 @@ subroutine gsd_limit_ocean_q(qinc,it)
   deallocate(rhgues)
 end subroutine gsd_limit_ocean_q 
 
-subroutine gsd_update_th2(tinc,it)
+subroutine gsd_update_t2m(tinc,it)
 !$$$  subprogram documentation block
 !                .      .    .                                       .
-! subprogram:    gsd_update_th2    adjust 2-m t based on analysis increment
+! subprogram:    gsd_update_t2m    adjust 2-m t based on analysis increment
 !   prgmmr: Hu          org: GSD                date: 2011-10-04
 !
 ! abstract:  This routine does the following things:
@@ -516,7 +516,8 @@ subroutine gsd_update_th2(tinc,it)
   use jfunc, only:  tsensible
   use constants, only: zero,one,fv,rd_over_cp_mass,one_tenth
   use gridmod, only: lat2,lon2,aeta1_ll,pt_ll,aeta2_ll
-! use guess_grids, only: nfldsig
+  use guess_grids, only:ges_prsi
+  use constants, only: half
 
   implicit none
 
@@ -530,14 +531,14 @@ subroutine gsd_update_th2(tinc,it)
   real(r_kind) :: dth2, work_prsl,work_prslk
 
   real(r_kind),dimension(:,:  ),pointer:: ges_ps =>NULL()
-  real(r_kind),dimension(:,:  ),pointer:: ges_th2=>NULL()
+  real(r_kind),dimension(:,:  ),pointer:: ges_t2m=>NULL()
   real(r_kind),dimension(:,:,:),pointer:: ges_q  =>NULL()
 
 !*******************************************************************************
 !
 ! 2-m temperature
 !  do it=1,nfldsig
-     call gsi_bundlegetpointer(gsi_metguess_bundle(it),'th2m',ges_th2,ier)
+     call gsi_bundlegetpointer(gsi_metguess_bundle(it),'t2m',ges_t2m,ier)
      if(ier/=0) return
      call gsi_bundlegetpointer(gsi_metguess_bundle(it),'ps',ges_ps,ier)
      if(ier/=0) return
@@ -552,17 +553,15 @@ subroutine gsd_update_th2(tinc,it)
               if(ihaveq/=0) cycle
               dth2=tinc(i,j)/(one+fv*ges_q(i,j,1))
            endif
-!          Convert sensible temperature to potential temperature
-           work_prsl  = one_tenth*(aeta1_ll(1)*(r10*ges_ps(i,j)-pt_ll)+ &
-                                   aeta2_ll(1) + pt_ll)
-           work_prslk = (work_prsl/r100)**rd_over_cp_mass
-           ges_th2(i,j) = ges_th2(i,j) + dth2/work_prslk
+
+!        do not need to convert sensible temperature to potential temperature
+           ges_t2m(i,j) = ges_t2m(i,j) + dth2
         end do
      end do
 !  end do
 
   return
-end subroutine gsd_update_th2
+end subroutine gsd_update_t2m
 
 subroutine gsd_update_q2(qinc,it)
 !$$$  subprogram documentation block

src/gsi/setupt.f90

@@ -1,5 +1,4 @@
 module t_setup
-
   implicit none
   private
   public:: setup
@@ -30,7 +29,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
   use m_obsdiagNode, only: obsdiagNode_assert
 
   use obsmod, only: sfcmodel,perturb_obs,oberror_tune,lobsdiag_forenkf,ianldate,&
-       lobsdiagsave,nobskeep,lobsdiag_allocated,time_offset,aircraft_recon 
+       lobsdiagsave,nobskeep,lobsdiag_allocated,time_offset,aircraft_recon
   use m_obsNode, only: obsNode
   use m_tNode, only: tNode
   use m_tNode, only: tNode_appendto
@@ -70,7 +69,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
   use converr, only: ptabl
   use rapidrefresh_cldsurf_mod, only: l_gsd_terrain_match_surftobs,l_sfcobserror_ramp_t
   use rapidrefresh_cldsurf_mod, only: l_pbl_pseudo_surfobst, pblh_ration,pps_press_incr
-  use rapidrefresh_cldsurf_mod, only: i_use_2mt4b,i_sfct_gross,l_closeobs,i_coastline       
+  use rapidrefresh_cldsurf_mod, only: i_use_2mt4b,i_sfct_gross,l_closeobs,i_coastline    
 
   use aircraftinfo, only: npredt,predt,aircraft_t_bc_pof,aircraft_t_bc, &
        aircraft_t_bc_ext,ostats_t,rstats_t,upd_pred_t
@@ -316,7 +315,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
   type(obs_diags),pointer:: my_diagLL
 
   real(r_kind) :: thisPBL_height,ratio_PBL_height,prestsfc,diffsfc,dthetav
-  real(r_kind) :: tges2m,qges2m,tges2m_water,qges2m_water, zges
+  real(r_kind) :: tges2m,qges2m,tges2m_water,qges2m_water
   real(r_kind) :: hr_offset
 
   equivalence(rstation_id,station_id)
@@ -329,8 +328,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
   real(r_kind),allocatable,dimension(:,:,:,:) :: ges_tv
   real(r_kind),allocatable,dimension(:,:,:,:) :: ges_q
   real(r_kind),allocatable,dimension(:,:,:  ) :: ges_q2
-  real(r_kind),allocatable,dimension(:,:,:  ) :: ges_t2
-  real(r_kind),allocatable,dimension(:,:,:  ) :: ges_th2
+  real(r_kind),allocatable,dimension(:,:,:  ) :: ges_t2m
 
   logical:: l_pbl_pseudo_itype
   integer(i_kind):: ich0
@@ -365,7 +363,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
 
 !  call GSD terrain match for surface temperature observation
   if(l_gsd_terrain_match_surftobs) then
-     call gsd_terrain_match_surfTobs(mype,nele,nobs,data) 
+     call gsd_terrain_match_surfTobs(mype,nele,nobs,data)
   endif
 
 !    index information for data array (see reading routine)
@@ -674,27 +672,18 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
      end if
 
 ! Interpolate log(ps) & log(pres) at mid-layers to obs locations/times
-     ! for global_2mDA case,  skipping reading the pressure levels into the met guess (to speed up the observer) 
-     ! these are not currently used here, so skip for now.
-
-     if (.not. do_global_2mDA) then  
-         call tintrp2a11(ges_ps,psges,dlat,dlon,dtime,hrdifsig,&
-              mype,nfldsig)
-         call tintrp2a1(ges_lnprsl,prsltmp,dlat,dlon,dtime,hrdifsig,&
-              nsig,mype,nfldsig)
-
-         drpx=zero
-         !if(sfctype .and. .not.twodvar_regional .and. .not. do_global_2mDA) then
-         if(sfctype .and. .not.twodvar_regional ) then
-            drpx=abs(one-((one/exp(dpres-log(psges))))**rd_over_cp)*t0c
-         end if
-
-    !    Put obs pressure in correct units to get grid coord. number
-         call grdcrd1(dpres,prsltmp(1),nsig,-1)
-     else 
-        drpx = zero 
-        dpres = one  ! put obs at surface
-     endif 
+     call tintrp2a11(ges_ps,psges,dlat,dlon,dtime,hrdifsig,&
+          mype,nfldsig)
+     call tintrp2a1(ges_lnprsl,prsltmp,dlat,dlon,dtime,hrdifsig,&
+          nsig,mype,nfldsig)
+
+     drpx=zero
+     if(sfctype .and. .not.twodvar_regional) then
+        drpx=abs(one-((one/exp(dpres-log(psges))))**rd_over_cp)*t0c
+     end if
+
+!    Put obs pressure in correct units to get grid coord. number
+     call grdcrd1(dpres,prsltmp(1),nsig,-1)
 
 ! Implementation of forward model ----------
 
@@ -732,14 +721,14 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
      elseif (sfctype .and. do_global_2mDA) then 
 ! SCENARIO 2: obs is sfctype, and do_global_2mDA scheme is on. 
 ! 2m forecast has been read from the sfc guess files
-! note: any changes to ges_t2 and tob in this block should also be made in buddy_check_t 
+! note: any changes to ges_t2m and tob in this block should also be made in buddy_check_t 
 
     ! mask: 0 - sea, 1 - land, 2-ice (val + 3 = interpolated from at least one different land cover
     ! for now, use only pure land (interpolation from mixed grid cells is pretty sketchy, but 
     ! may result in lack of obs near coasts) 
           if (int(data(idomsfc,i)) .NE. 1  ) muse(i) = .false. 
 
-          call tintrp2a11(ges_t2,tges2m,dlat,dlon,dtime,hrdifsig,&
+          call tintrp2a11(ges_t2m,tges2m,dlat,dlon,dtime,hrdifsig,&
             mype,nfldsig)
 
 ! correct obs to model terrain height - equivalent to gsd_terrain_match, which  uses lapse rate 
@@ -752,7 +741,16 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
           !update the station elevation
           data(istnelv,i) = data(izz,i) 
 
-
+          if(save_jacobian) then
+             t2m_ind = getindex(svars2d, 't2m')
+             if (t2m_ind < 0) then
+                 print *, 'Error: no variable t2m in state vector.Exiting.'
+                 call stop2(1300)
+             endif
+             dhx_dx%st_ind(1) = sum(levels(1:ns3d))  + t2m_ind
+             dhx_dx%end_ind(1) = sum(levels(1:ns3d)) + t2m_ind
+             dhx_dx%val(1) = one
+          endif
      else
 ! SCENARIO 3: obs is sfctype, and neither sfcmodel nor do_global_2mDA is chosen
 !        .or. obs is not sfctype     
@@ -805,13 +803,12 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
 
 ! SCENARIO 4: obs is sfctype, and i_use_2mt4b flag is on (turns on LAM sfc DA)
         if(i_use_2mt4b>0 .and. sfctype) then
+
            if(i_coastline==1 .or. i_coastline==3) then
 
 !          Interpolate guess th 2m to observation location and time
-              call tintrp2a11_csln(ges_th2,tges2m,tges2m_water,dlat,dlon,dtime,hrdifsig,&
+              call tintrp2a11_csln(ges_t2m,tges2m,tges2m_water,dlat,dlon,dtime,hrdifsig,&
                 mype,nfldsig)
-              tges2m=tges2m*(r10*psges/r1000)**rd_over_cp_mass  ! convert to sensible T         
-              tges2m_water=tges2m_water*(r10*psges/r1000)**rd_over_cp_mass  ! convert to sensible T         
               if(iqtflg)then
                  call tintrp2a11_csln(ges_q2,qges2m,qges2m_water,dlat,dlon,dtime,hrdifsig,&
                      mype,nfldsig)
@@ -821,9 +818,8 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
               if( abs(tob-tges2m) > abs(tob-tges2m_water)) tges2m=tges2m_water
            else
 !          Interpolate guess th 2m to observation location and time
-              call tintrp2a11(ges_th2,tges2m,dlat,dlon,dtime,hrdifsig,&
+              call tintrp2a11(ges_t2m,tges2m,dlat,dlon,dtime,hrdifsig,&
                 mype,nfldsig)
-              tges2m=tges2m*(r10*psges/r1000)**rd_over_cp_mass  ! convert to sensible T         
               if(iqtflg)then
                  call tintrp2a11(ges_q2,qges2m,dlat,dlon,dtime,hrdifsig,&
                      mype,nfldsig)
@@ -856,6 +852,8 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
          endif
      else
         rlow  = zero
+        ramp  = zero
+        drpx  = zero
      endif
 
      rhgh=max(zero,dpres-rsigp-r0_001)
@@ -888,10 +886,8 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
 
 
 ! Compute innovation
-     !if(i_use_2mt4b>0 .and. sfctype) then
      if(sfctype .and. (( i_use_2mt4b>0) .or. do_global_2mDA) ) then
         ddiff = tob-tges2m
-        tges = tges2m ! not necessary? 
      else
         ddiff = tob-tges
      endif
@@ -1498,47 +1494,25 @@ subroutine init_vars_
          write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus
          call stop2(999)
      endif
-     if(i_use_2mt4b>0) then
-!    get th2m ...
-        varname='th2m'
-        call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus)
-        if (istatus==0) then
-            if(allocated(ges_th2))then
-               write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc '
-               call stop2(999)
-            endif
-            allocate(ges_th2(size(rank2,1),size(rank2,2),nfldsig))
-            ges_th2(:,:,1)=rank2
-            do ifld=2,nfldsig
-               call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank2,istatus)
-               ges_th2(:,:,ifld)=rank2
-            enddo
-        else
-            write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus
-            call stop2(999)
-        endif
-     endif
-     if( do_global_2mDA) then
+     if(i_use_2mt4b>0 .or. do_global_2mDA) then
 !    get t2m ...
         varname='t2m'
         call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus)
         if (istatus==0) then
-            if(allocated(ges_t2))then
+            if(allocated(ges_t2m))then
                write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc '
                call stop2(999)
             endif
-            allocate(ges_t2(size(rank2,1),size(rank2,2),nfldsig))
-            ges_t2(:,:,1)=rank2
+            allocate(ges_t2m(size(rank2,1),size(rank2,2),nfldsig))
+            ges_t2m(:,:,1)=rank2
             do ifld=2,nfldsig
                call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank2,istatus)
-               ges_t2(:,:,ifld)=rank2
+               ges_t2m(:,:,ifld)=rank2
             enddo
         else
             write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus
             call stop2(999)
         endif
-     endif
-     if(i_use_2mt4b>0 .or. do_global_2mDA) then
 !    get q2m ...
         varname='q2m'
         call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus)
@@ -1888,8 +1862,7 @@ subroutine final_vars_
     if(allocated(ges_u )) deallocate(ges_u )
     if(allocated(ges_ps)) deallocate(ges_ps)
     if(allocated(ges_q2)) deallocate(ges_q2)
-    if(allocated(ges_th2)) deallocate(ges_th2)
-    if(allocated(ges_t2)) deallocate(ges_t2)
+    if(allocated(ges_t2m)) deallocate(ges_t2m)
   end subroutine final_vars_
 
 end subroutine setupt

src/gsi/update_guess.f90

@@ -113,7 +113,7 @@ subroutine update_guess(sval,sbias)
   use mpimod, only: mype
   use constants, only: zero,one,fv,max_varname_length,qmin,qcmin,tgmin,&
                        r100,one_tenth,tiny_r_kind
-  use jfunc, only: iout_iter,bcoption,tsensible,clip_supersaturation,superfact
+  use jfunc, only: iout_iter,bcoption,tsensible,clip_supersaturation,superfact,do_global_2mDA
   use gridmod, only: lat2,lon2,nsig,&
        regional,twodvar_regional,regional_ozone,&
        l_reg_update_hydro_delz
@@ -136,7 +136,7 @@ subroutine update_guess(sval,sbias)
   use rapidrefresh_cldsurf_mod, only: l_gsd_limit_ocean_q,l_gsd_soilTQ_nudge
   use rapidrefresh_cldsurf_mod, only: i_use_2mq4b,i_use_2mt4b
   use gsd_update_mod, only: gsd_limit_ocean_q,gsd_update_soil_tq,&
-       gsd_update_th2,gsd_update_q2
+       gsd_update_t2m,gsd_update_q2
   use qcmod, only: pvis,pcldch,vis_thres,cldch_thres 
   use obsmod, only: l_wcp_cwm
   use directDA_radaruse_mod, only: l_use_cvpqx, l_use_dbz_directDA
@@ -454,15 +454,15 @@ subroutine update_guess(sval,sbias)
         endif
         call  gsd_update_soil_tq(tinc_1st,is_t,qinc_1st,is_q,it)
      endif  ! l_gsd_soilTQ_nudge
-     if (i_use_2mt4b > 0 .and. is_t>0) then
+     if ( (i_use_2mt4b > 0.or.do_global_2mDA) .and. is_t>0) then
         do j=1,lon2
            do i=1,lat2
               tinc_1st(i,j)=p_tv(i,j,1)
            end do
         end do
-        call  gsd_update_th2(tinc_1st,it)
-     endif ! l_gsd_th2_adjust
-     if (i_use_2mq4b > 0 .and. is_q>0) then
+        call  gsd_update_t2m(tinc_1st,it)
+     endif ! l_gsd_t2m_adjust
+     if ( (i_use_2mq4b > 0.or.do_global_2mDA) .and. is_q>0) then
         do j=1,lon2
            do i=1,lat2
               qinc_1st(i,j)=p_q(i,j,1)