Changeset 899 in lmdz_wrf

Timestamp:

Jun 19, 2016, 2:48:04 PM (9 years ago)

Author:

lfita

Message:

Adding `p_interp.F90' vertical presssure interpolation from C. Bruyere NCAR

File:

• trunk/tools/module_ForInterpolate.F90 (modified) (2 diffs)

trunk/tools/module_ForInterpolate.F90

@@ -1 @@
-! Module to interpolate values from a giving projection
+! Module to interpolate values from a giving projection and pressure interpolation
 ! To be included in a python
-! f2py -m module_ForInterpolate --f90exec=/usr/bin/gfortran-4.7 -c module_ForInterpolate.F90 module_generic.F90 >& run_f2py.log
+! f2py -m module_ForInterpolate --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForInterpolate.F90  >& run_f2py.log
 MODULE module_ForInterpolate
 
@@ -922 +922 @@
 END SUBROUTINE Interpolate1DLl
 
+
+ SUBROUTINE interp (data_out, data_in, pres_field, interp_levels, psfc, ter, tk, qv, ix, iy, iz, it, &
+                     num_metgrid_levels, LINLOG, extrapolate, GEOPT, MISSING)
+! Interpolation subroutine from the p_interp.F90 NCAR program
+!  Program to read wrfout data and interpolate to pressure levels
+!  The program reads namelist.pinterp
+!  November 2007 - Cindy Bruyere
+!
+     INTEGER                                          :: ix, iy, iz, it
+     INTEGER                                          :: num_metgrid_levels, LINLOG
+     REAL, DIMENSION(ix, iy, num_metgrid_levels, it)  :: data_out
+     REAL, DIMENSION(ix, iy, iz, it)                  :: data_in, pres_field, tk, qv
+     REAL, DIMENSION(ix, iy, it)                      :: psfc
+     REAL, DIMENSION(ix, iy)                          :: ter
+     REAL, DIMENSION(num_metgrid_levels)              :: interp_levels
+
+     INTEGER                                          :: i, j, itt, k, kk, kin
+     REAL, DIMENSION(num_metgrid_levels)              :: data_out1D
+     REAL, DIMENSION(iz)                              :: data_in1D, pres_field1D
+     INTEGER                                          :: extrapolate
+     REAL                                             :: MISSING
+     REAL, DIMENSION(ix, iy, num_metgrid_levels, it)  :: N
+     REAL                                             :: sumA, sumN, AVE_geopt
+     LOGICAL                                          :: GEOPT
+
+!!!!!!! Variables
+! data_out: interpolated field
+! data_in: field to interpolate
+! pres_field: pressure field [Pa]
+! interp_levels: pressure levels to interpolate [hPa]
+! psfc: surface pressure [Pa]
+! ter: terrein height [m]
+! tk: temperature [K]
+! qv: mositure mizing ratio [kg/kg]
+! i[x/y/z/t]: size of the matrices
+! num_metgrid_levels: number of pressure values to interpolate
+! LINLOG: if abs(linlog)=1 use linear interp in pressure
+!         if abs(linlog)=2 linear interp in ln(pressure)
+! extrapolate: whether to set to missing value below/above model ground and top (0), or extrapolate (1)
+! GEOPT: Wether the file is the geopotential file or not
+! MISSING: Missing value
+
+     N = 1.0
+
+     expon=287.04*.0065/9.81
+
+     do itt = 1, it
+        do j = 1, iy
+        do i = 1, ix
+           data_in1D(:)    = data_in(i,j,:,itt)
+           pres_field1D(:) = pres_field(i,j,:,itt)
+           CALL int1D (data_out1D, data_in1D, pres_field1D, interp_levels, iz, num_metgrid_levels, LINLOG, MISSING)
+           data_out(i,j,:,itt) = data_out1D(:)
+        end do
+        end do
+     end do
+
+
+     ! Fill in missing values
+     IF ( extrapolate == 0 ) RETURN       !! no extrapolation - we are out of here
+
+     ! First find where about 400 hPa is located
+     kk = 0
+     find_kk : do k = 1, num_metgrid_levels
+        kk = k
+        if ( interp_levels(k) <= 40000. ) exit find_kk
+     end do find_kk
+
+     
+     IF ( GEOPT ) THEN     !! geopt is treated different below ground
+
+        do itt = 1, it
+           do k = 1, kk
+              do j = 1, iy
+              do i = 1, ix
+                 IF ( data_out(i,j,k,itt) == MISSING .AND. interp_levels(k) < psfc(i,j,itt) ) THEN
+
+!                We are below the first model level, but above the ground 
+
+                    data_out(i,j,k,itt) = ((interp_levels(k) - pres_field(i,j,1,itt))*ter(i,j)*9.81 +  &
+                                           (psfc(i,j,itt) - interp_levels(k))*data_in(i,j,1,itt) ) /   &
+                                          (psfc(i,j,itt) - pres_field(i,j,1,itt))
+
+                 ELSEIF ( data_out(i,j,k,itt) == MISSING ) THEN
+
+!                We are below both the ground and the lowest data level.
+
+!                First, find the model level that is closest to a "target" pressure
+!                level, where the "target" pressure is delta-p less that the local
+!                value of a horizontally smoothed surface pressure field.  We use
+!                delta-p = 150 hPa here. A standard lapse rate temperature profile
+!                passing through the temperature at this model level will be used
+!                to define the temperature profile below ground.  This is similar
+!                to the Benjamin and Miller (1990) method, except that for
+!                simplicity, they used 700 hPa everywhere for the "target" pressure.
+!                Code similar to what is implemented in RIP4
+
+                    ptarget = (psfc(i,j,itt)*.01) - 150.
+                    dpmin=1.e4
+                    kupper = 0
+                    loop_kIN : do kin=iz,1,-1
+                       kupper = kin
+                       dp=abs( (pres_field(i,j,kin,itt)*.01) - ptarget )
+                       if (dp.gt.dpmin) exit loop_kIN
+                       dpmin=min(dpmin,dp)
+                    enddo loop_kIN
+
+                    pbot=max(pres_field(i,j,1,itt),psfc(i,j,itt))
+                    zbot=min(data_in(i,j,1,itt)/9.81,ter(i,j))
+
+                    tbotextrap=tk(i,j,kupper,itt)*(pbot/pres_field(i,j,kupper,itt))**expon
+                    tvbotextrap=virtual(tbotextrap,qv(i,j,1,itt))
+
+                    data_out(i,j,k,itt) = (zbot+tvbotextrap/.0065*(1.-(interp_levels(k)/pbot)**expon))*9.81
+               
+                 ENDIF
+              enddo
+              enddo
+           enddo
+        enddo
+
+
+        !!! Code for filling missing data with an average - we don't want to do this
+        !!do itt = 1, it
+           !!loop_levels : do k = 1, num_metgrid_levels
+              !!sumA = SUM(data_out(:,:,k,itt), MASK = data_out(:,:,k,itt) /= MISSING)
+              !!sumN = SUM(N(:,:,k,itt), MASK = data_out(:,:,k,itt) /= MISSING)
+              !!IF ( sumN == 0. ) CYCLE loop_levels
+              !!AVE_geopt = sumA/sumN
+              !!WHERE ( data_out(:,:,k,itt) == MISSING )
+                 !!data_out(:,:,k,itt) = AVE_geopt
+              !!END WHERE
+           !!end do loop_levels
+        !!end do
+
+     END IF
+     
+     !!! All other fields and geopt at higher levels come here
+     do itt = 1, it
+        do j = 1, iy
+        do i = 1, ix
+          do k = 1, kk
+             if ( data_out(i,j,k,itt) == MISSING ) data_out(i,j,k,itt) = data_in(i,j,1,itt)
+          end do
+          do k = kk+1, num_metgrid_levels
+             if ( data_out(i,j,k,itt) == MISSING ) data_out(i,j,k,itt) = data_in(i,j,iz,itt)
+          end do
+        end do
+        end do
+     end do
+
+ END SUBROUTINE interp 
+
+ SUBROUTINE int1D(xxout, xxin, ppin, ppout, npin, npout, LINLOG, MISSING)
+
+! Modified from int2p - NCL code
+! routine to interpolate from one set of pressure levels
+! .   to another set  using linear or ln(p) interpolation
+!
+! NCL: xout = int2p (pin,xin,pout,linlog)
+! This code was originally written for a specific purpose.
+! .   Several features were added for incorporation into NCL's
+! .   function suite including linear extrapolation.
+!
+! nomenclature:
+!
+! .   ppin   - input pressure levels. The pin can be
+! .            be in ascending or descending order
+! .   xxin   - data at corresponding input pressure levels
+! .   npin   - number of input pressure levels >= 2
+! .   ppout  - output pressure levels (input by user)
+! .            same (ascending or descending) order as pin
+! .   xxout  - data at corresponding output pressure levels
+! .   npout  - number of output pressure levels
+! .   linlog - if abs(linlog)=1 use linear interp in pressure
+! .            if abs(linlog)=2 linear interp in ln(pressure)
+! .   missing- missing data code. 
+
+!                                                ! input types
+      INTEGER   :: npin,npout,linlog,ier
+      real      :: ppin(npin),xxin(npin),ppout(npout)
+      real      :: MISSING       
+     logical                                          :: AVERAGE
+!                                                ! output
+      real      :: xxout(npout)
+      INTEGER   :: j1,np,nl,nin,nlmax,nplvl
+      INTEGER   :: nlsave,np1,no1,n1,n2,nlstrt
+      real      :: slope,pa,pb,pc
+
+! automatic arrays
+      real      :: pin(npin),xin(npin),p(npin),x(npin)
+      real      :: pout(npout),xout(npout)
+
+
+      xxout = MISSING
+      pout  = ppout
+      p     = ppin
+      x     = xxin
+      nlmax = npin
+
+! exact p-level matches
+      nlstrt = 1
+      nlsave = 1
+      do np = 1,npout
+          xout(np) = MISSING
+          do nl = nlstrt,nlmax
+              if (pout(np).eq.p(nl)) then
+                  xout(np) = x(nl)
+                  nlsave = nl + 1
+                  go to 10
+              end if
+          end do
+   10     nlstrt = nlsave
+      end do
+
+      if (LINLOG.eq.1) then
+          do np = 1,npout
+              do nl = 1,nlmax - 1
+                  if (pout(np).lt.p(nl) .and. pout(np).gt.p(nl+1)) then
+                      slope = (x(nl)-x(nl+1))/ (p(nl)-p(nl+1))
+                      xout(np) = x(nl+1) + slope* (pout(np)-p(nl+1))
+                  end if
+              end do
+          end do
+      elseif (LINLOG.eq.2) then
+          do np = 1,npout
+              do nl = 1,nlmax - 1
+                  if (pout(np).lt.p(nl) .and. pout(np).gt.p(nl+1)) then
+                      pa = log(p(nl))
+                      pb = log(pout(np))
+! special case: in case someone inadvertently enter p=0.
+                      if (p(nl+1).gt.0.d0) then
+                          pc = log(p(nl+1))
+                      else
+                          pc = log(1.d-4)
+                      end if
+
+                      slope = (x(nl)-x(nl+1))/ (pa-pc)
+                      xout(np) = x(nl+1) + slope* (pb-pc)
+                  end if
+              end do
+          end do
+      end if
+
+
+! place results in the return array;
+      xxout = xout
+
+ END SUBROUTINE int1D
+
 END MODULE module_ForInterpolate