Added the capability for dynamic lapse rate adjustment

This capability requires each forcing reader to specify the dynamic
lapse rates to be used. Currently implemented for MERRA2

Resolves NASA-LIS#410

lis/core/LIS_metforcingMod.F90

@@ -95,6 +95,7 @@ module LIS_metforcingMod
 
   type, public :: forc_dec_type
      real, allocatable :: modelelev(:)
+     real, allocatable :: lapseRate(:)
   end type forc_dec_type
 
 !BOP
@@ -168,6 +169,10 @@ subroutine metforcing_init_offline
                 write(LIS_logunit,*) "[INFO] Elevation correction turned on for:  ",&
                 trim(LIS_rc%metforc(m))
                 allocate(LIS_forc(n,m)%modelelev(LIS_rc%ngrid(n)))
+                allocate(LIS_forc(n,m)%lapseRate(LIS_rc%ngrid(n)))
+
+                !initialized to the uniform value
+                LIS_forc(n,m)%lapseRate = -0.0065
              endif
           enddo
        enddo
@@ -1101,7 +1106,9 @@ subroutine LIS_get_met_forcing(n)
                 call LIS_endrun
              endif
 
-             call LIS_lapseRateCorrection(n,LIS_forc(n,m)%modelelev,&
+             call LIS_lapseRateCorrection(n,&
+                  LIS_forc(n,m)%modelelev,&
+                  LIS_forc(n,m)%lapseRate,&
                   LIS_FORC_Base_State(n,m))
 
           elseif(LIS_rc%met_ecor(m).eq."slope-aspect") then 
@@ -1129,7 +1136,9 @@ subroutine LIS_get_met_forcing(n)
                 call LIS_endrun
              endif
 
-             call LIS_lapseRateCorrection(n, LIS_forc(n,m)%modelelev,&
+             call LIS_lapseRateCorrection(n, &
+                  LIS_forc(n,m)%modelelev,&
+                  LIS_forc(n,m)%lapseRate,&
                   LIS_FORC_Base_State(n,m))
              call LIS_slopeAspectCorrection(n, LIS_FORC_Base_State(n,m))
 

lis/core/LIS_spatialDownscalingMod.F90

@@ -71,13 +71,14 @@ module LIS_spatialDownscalingMod
 !  14 Nov 2003: Sujay Kumar; Adopted in LIS
 !
 ! !INTERFACE:
-  subroutine LIS_lapseRateCorrection(nest, modelelev, LIS_FORC_Base_State)
+  subroutine LIS_lapseRateCorrection(nest, modelelev, lapseRate, LIS_FORC_Base_State)
 ! !USES:
 
     implicit none
 ! !ARGUMENTS: 
     integer, intent(in) :: nest
     real                :: modelelev(LIS_rc%ngrid(nest))
+    real                :: lapseRate(LIS_rc%ngrid(nest))
     type(ESMF_State)    :: LIS_FORC_Base_State
 
 ! !DESCRIPTION:
@@ -116,7 +117,7 @@ subroutine LIS_lapseRateCorrection(nest, modelelev, LIS_FORC_Base_State)
     real, pointer      :: tmp(:),q2(:),lwd(:),psurf(:)
 
     rdry = 287.
-    lapse = -0.0065
+!    lapse = -0.0065
 
     call ESMF_StateGet(LIS_FORC_Base_State,&
          trim(LIS_FORC_Tair%varname(1)),tmpField,&
@@ -155,12 +156,13 @@ subroutine LIS_lapseRateCorrection(nest, modelelev, LIS_FORC_Base_State)
 
 
     do t=1,LIS_rc%ntiles(nest)
-       if(tmp(t).gt.0) then 
+       if(tmp(t).gt.0) then
           force_tmp = tmp(t)
           force_hum = q2(t)
           force_lwd = lwd(t)
           force_prs = psurf(t)
           index = LIS_domain(nest)%tile(t)%index
+          lapse = lapseRate(index)
           elevdiff = LIS_domain(nest)%tile(t)%elev-&
                modelelev(index)
           tcforce=force_tmp+(lapse*elevdiff)

lis/metforcing/merra2/get_merra2.F90

@@ -174,6 +174,9 @@ subroutine get_merra2(n, findex)
   integer           :: hr_int1, hr_int2
   integer           :: movetime  ! Flag to move bookend2 files to bookend1
 
+  character(len=LIS_CONST_PATH_LEN) :: lapseratefname
+  character*8                       :: fdate
+
 ! _________________________________________________________
 
 ! Please note that the timing logic has been tested only for
@@ -307,8 +310,14 @@ subroutine get_merra2(n, findex)
      do kk= merra2_struc(n)%st_iterid, merra2_struc(n)%en_iterid
         call merra2files(n,kk,findex,merra2_struc(n)%merra2dir, yr1, mo1, da1, &
              slvname, flxname, lfoname, radname)
+
+        write(unit=fdate, fmt='(i4.4,i2.2,i2.2)') yr1,mo1,da1
+        lapseratefname = trim(merra2_struc(n)%dynlapseratedir)//&
+             '/MERRA2.lapse_rate.hourly.'//&
+             trim(fdate)//'.global.nc'
+
         call read_merra2(n, order, mo1, &
-             findex, slvname, flxname, lfoname, radname,&
+             findex, slvname, flxname, lfoname, radname, lapseratefname,&
              merra2_struc(n)%merraforc1(kk,:,:,:), ferror)
      enddo
   endif
@@ -325,8 +334,15 @@ subroutine get_merra2(n, findex)
      do kk= merra2_struc(n)%st_iterid, merra2_struc(n)%en_iterid
         call merra2files(n,kk,findex,merra2_struc(n)%merra2dir, yr2, mo2, da2, &
              slvname, flxname, lfoname, radname)
+
+        write(unit=fdate, fmt='(i4.4,i2.2,i2.2)') yr2,mo2,da2
+        lapseratefname = trim(merra2_struc(n)%dynlapseratedir)//&
+             '/MERRA2.lapse_rate.hourly.'//&
+             trim(fdate)//'.global.nc'
+
+
         call read_merra2(n, order, mo2,&
-             findex, slvname, flxname, lfoname, radname, &
+             findex, slvname, flxname, lfoname, radname, lapseratefname,&
              merra2_struc(n)%merraforc2(kk,:,:,:), ferror)
      enddo
   endif

lis/metforcing/merra2/merra2_forcingMod.F90

@@ -125,6 +125,10 @@ module merra2_forcingMod
      real, allocatable       :: merraxrange(:,:,:,:)
      real, allocatable       :: merracdf(:,:,:,:)
      integer, allocatable    :: rseed(:,:)
+
+     integer                 :: usedynlapserate
+     character(len=LIS_CONST_PATH_LEN) :: dynlapseratedir
+
   end type merra2_type_dec
 
   type(merra2_type_dec), allocatable :: merra2_struc(:)

lis/metforcing/merra2/read_merra2.F90

@@ -17,9 +17,10 @@
 ! 18 Mar 2015: James Geiger, initial code (based on merra-land)
 !
 ! !INTERFACE:
-subroutine read_merra2(n, order, month, findex,          &
-                       slvname, flxname, lfoname, radname, &
-                       merraforc, ferror)
+subroutine read_merra2(n, order, month, findex,  &
+     slvname, flxname, lfoname, radname,         &
+     lapseratefname,                             &
+     merraforc, ferror)
 ! !USES:
   use LIS_coreMod,       only : LIS_rc, LIS_domain, LIS_masterproc
   use LIS_logMod
@@ -40,6 +41,7 @@ subroutine read_merra2(n, order, month, findex,          &
   character(len=*), intent(in) :: flxname
   character(len=*), intent(in) :: lfoname
   character(len=*), intent(in) :: radname
+  character(len=*), intent(in) :: lapseratefname
   real, intent(inout)          :: merraforc(merra2_struc(n)%nvars, 24, &
        LIS_rc%lnc(n)*LIS_rc%lnr(n))
   integer, intent(out)         :: ferror          
@@ -113,7 +115,15 @@ subroutine read_merra2(n, order, month, findex,          &
   real      :: pardr(merra2_struc(n)%ncold, merra2_struc(n)%nrold,24)
   real      :: pardf(merra2_struc(n)%ncold, merra2_struc(n)%nrold,24)
   real      :: hlml(merra2_struc(n)%ncold, merra2_struc(n)%nrold,24)
-!  real      :: emis(merra2_struc(n)%ncold, merra2_struc(n)%nrold,24)
+  !  real      :: emis(merra2_struc(n)%ncold, merra2_struc(n)%nrold,24)
+
+
+  integer                           :: ftn_drate
+  integer                           :: gid
+  integer                           :: lapserateid
+  character*8                       :: fdate
+  real, allocatable                 :: lapse_rate_in(:,:)
+  real, allocatable                 :: lapse_rate_out(:)
 ! __________________________________________________________________________
 
 #if (defined USE_NETCDF3) 
@@ -214,6 +224,54 @@ subroutine read_merra2(n, order, month, findex,          &
      call interp_merra2_var(n,findex,month,uwind, 5, .false., merraforc)
      call interp_merra2_var(n,findex,month,vwind, 6, .false., merraforc)
      call interp_merra2_var(n,findex,month,ps,    7, .false., merraforc)
+
+     if(merra2_struc(n)%usedynlapserate.eq.1) then
+
+        inquire(file=trim(lapseratefname),exist=file_exists)
+
+        if(file_exists) then
+           allocate(lapse_rate_in(merra2_struc(n)%ncold, merra2_struc(n)%nrold))
+           allocate(lapse_rate_out(LIS_rc%lnc(n)*LIS_rc%lnr(n)))
+
+           call LIS_verify(nf90_open(path=trim(lapseratefname),&
+                mode=NF90_NOWRITE,&
+                ncid = ftn_drate),&
+                'nf90_open failed for '//trim(lapseratefname))
+
+           call LIS_verify(nf90_inq_varid(ftn_drate,&
+                'lapse_rate',lapserateid),&
+                'nf90_inq_varid failed for lapse_rate')
+           call LIS_verify(nf90_get_var(ftn_drate,&
+                lapserateid,lapse_rate_in),&
+                'nf90_get_var failed for lapse_rate')
+           call LIS_verify(nf90_close(ftn_drate),&
+                'nf90_close failed')
+
+
+           call interp_lapserate_var(n,&
+                lapse_rate_in,&
+                lapse_rate_out)                
+
+           do r=1,LIS_rc%lnr(n)
+              do c=1,LIS_rc%lnc(n)
+                 if(LIS_domain(n)%gindex(c,r).ne.-1) then
+
+                    gid = LIS_domain(n)%gindex(c,r)                    
+                    LIS_forc(n,findex)%lapseRate(gid) = &
+                         lapse_rate_out(c+(r-1)*LIS_rc%lnc(n))/1000.0
+                 endif
+              enddo
+           enddo
+
+           deallocate(lapse_rate_in)
+           deallocate(lapse_rate_out)
+
+        else
+           LIS_forc(n,findex)%lapseRate(:) = -0.0065
+        endif
+
+     endif
+
   else
      write(LIS_logunit,*) '[ERR] ',trim(slvname)//' does not exist'
      call LIS_endrun()
@@ -717,6 +775,73 @@ subroutine interp_merra2_var(n,findex,month, input_var,  var_index, &
 
 end subroutine interp_merra2_var
 
+
+!BOP
+! 
+! !ROUTINE: interp_lapserate_var
+! \label{interp_lapserate_var}
+! 
+! !INTERFACE: 
+subroutine interp_lapserate_var(n,input_var, output_var)
+
+! !USES: 
+  use LIS_coreMod
+  use LIS_logMod
+  use LIS_spatialDownscalingMod
+  use merra2_forcingMod, only : merra2_struc
+#if(defined USE_NETCDF3 || defined USE_NETCDF4)      
+  use netcdf
+#endif
+  implicit none
+
+! !ARGUMENTS:
+  integer, intent(in)    :: n
+  real,    intent(in)    :: input_var(merra2_struc(n)%ncold, &
+       merra2_struc(n)%nrold)
+  real,    intent(inout) :: output_var(LIS_rc%lnc(n)*LIS_rc%lnr(n))
+  !
+! !DESCRIPTION: 
+!  This subroutine spatially interpolates a MERRA2 field
+!  to the LIS running domain
+! 
+!EOP
+
+  integer   :: t,c,r,k,iret
+  integer   :: doy
+  integer   :: ftn
+  real      :: f (merra2_struc(n)%ncold*merra2_struc(n)%nrold)
+  logical*1 :: lb(merra2_struc(n)%ncold*merra2_struc(n)%nrold)
+  logical*1 :: lo(LIS_rc%lnc(n)*LIS_rc%lnr(n))
+  integer   :: input_size
+
+! _____________________________________________________________
+
+  input_size = merra2_struc(n)%ncold*merra2_struc(n)%nrold
+
+  lb = .true.
+  do r=1,merra2_struc(n)%nrold
+     do c=1,merra2_struc(n)%ncold
+        k= c+(r-1)*merra2_struc(n)%ncold
+        f(k) = input_var(c,r)
+        if ( f(k) == 1.e+15 ) then 
+           f(k)  = LIS_rc%udef
+           lb(k) = .false. 
+        endif
+     enddo
+  enddo
+
+  call bilinear_interp(LIS_rc%gridDesc(n,:),lb,f,lo,&
+       output_var, &
+       merra2_struc(n)%mi,LIS_rc%lnc(n)*LIS_rc%lnr(n), & 
+       LIS_domain(n)%lat, LIS_domain(n)%lon,&
+       merra2_struc(n)%w111,merra2_struc(n)%w121,&
+       merra2_struc(n)%w211,merra2_struc(n)%w221,&
+       merra2_struc(n)%n111,merra2_struc(n)%n121,&
+       merra2_struc(n)%n211,merra2_struc(n)%n221,&
+       LIS_rc%udef, iret)
+
+end subroutine interp_lapserate_var
+
 #if 0 
 !BOP
 !

lis/metforcing/merra2/readcrd_merra2.F90

@@ -32,6 +32,7 @@ subroutine readcrd_merra2()
   implicit none
 
   integer :: n,t,rc
+  logical :: usedynlapserate
 
   call ESMF_ConfigFindLabel(LIS_config,"MERRA2 forcing directory:",rc=rc)
   do n=1,LIS_rc%nnest
@@ -59,6 +60,34 @@ subroutine readcrd_merra2()
           'MERRA2 use corrected total precipitation: not defined')
   enddo
 
+  call ESMF_ConfigFindLabel(LIS_config,"MERRA2 apply dynamic lapse rates:",rc=rc)
+  do n=1,LIS_rc%nnest
+     call ESMF_ConfigGetAttribute(LIS_config,merra2_struc(n)%usedynlapserate,&
+          default=0, rc=rc)
+     call LIS_verify(rc,&
+          'MERRA2 apply dynamic lapse rates: not defined')
+  enddo
+
+  usedynlapserate = .true.
+
+  do n=1,LIS_rc%nnest
+     if(merra2_struc(n)%usedynlapserate.eq.0) then
+        usedynlapserate = .false. 
+     endif
+  enddo
+
+  if(usedynlapserate) then
+     call ESMF_ConfigFindLabel(LIS_config,"MERRA2 dynamic lapse rate data directory:",rc=rc)
+  do n=1,LIS_rc%nnest
+     call ESMF_ConfigGetAttribute(LIS_config,merra2_struc(n)%dynlapseratedir,&
+          rc=rc)
+     call LIS_verify(rc,&
+          'MERRA2 dynamic lapse rate data directory: not defined')
+  enddo
+
+
+  endif
+
   do n=1,LIS_rc%nnest
      merra2_struc(n)%usescalef = 0 
      merra2_struc(n)%nbins = 50 !currently hardcoded