nilyin
/
wrf-python
forked from 3rdparty/wrf-python
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Removed ncl_reference directory since it is not needed anymore

main
Bill Ladwig 8 years ago
parent
17cfd4fbbd
commit
88b8623ac8
28 changed files with 0 additions and 33110 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 5058
      ncl_reference/WRFUserARW.ncl
  2. 293
      ncl_reference/WRF_contributed.ncl
  3. 158
      ncl_reference/calc_uh.f90
  4. 60
      ncl_reference/cloud_fracf.f
  5. 72
      ncl_reference/eqthecalc.f
  6. 82
      ncl_reference/etaconv.py
  7. 183
      ncl_reference/jinja_context.py
  8. 170
      ncl_reference/module_model_constants.f
  9. 77
      ncl_reference/product_table.txt
  10. 4575
      ncl_reference/psadilookup.dat
  11. 215
      ncl_reference/rcm2points.f
  12. 376
      ncl_reference/rcm2rgrid.f
  13. 832
      ncl_reference/rcmW.c
  14. 612
      ncl_reference/rip_cape.f
  15. 13253
      ncl_reference/wrfW.c
  16. 402
      ncl_reference/wrf_bint3d.f
  17. 763
      ncl_reference/wrf_cloud_topW.c
  18. 117
      ncl_reference/wrf_fctt.f
  19. 1917
      ncl_reference/wrf_fddaobs_in.F
  20. 380
      ncl_reference/wrf_map_fix.ncl
  21. 109
      ncl_reference/wrf_pvo.f
  22. 100
      ncl_reference/wrf_relhl.f
  23. 264
      ncl_reference/wrf_rip_phys_routines.f
  24. 771
      ncl_reference/wrf_user.f
  25. 209
      ncl_reference/wrf_user_dbz.f
  26. 511
      ncl_reference/wrf_user_latlon_routines.f
  27. 404
      ncl_reference/wrf_vinterp.f
  28. 1147
      ncl_reference/wrf_vinterpW.c

5058
ncl_reference/WRFUserARW.ncl
View File

File diff suppressed because it is too large Load Diff

293
ncl_reference/WRF_contributed.ncl
Unescape View File

@ -1,293 +0,0 @@
;*************************************************************************
; Note: several of the functions/procedures are used
; to invoke old [ugly!] function names.
;*************************************************************************
; D. Shea
;
; convert WRF character variable "Times" to udunits
; 2001-06-11_12:00:00
;
; convert WRF character variable "Times" to a coordinate variable "Time"
; opt can be "integer" or "string"
; . integer: opt = 0 : hours since initial time: Times(0,:)
; . opt = 1 : hours since 1901-01-01 00:00:00
; . string: opt = 'any udunits compatible string'
;
undef ("WRF_Times2Udunits_c")
function WRF_Times2Udunits_c(Times:character, opt)
local dimT, rank, year, month, day, hour, minute, sec, units, time
begin
dimT = dimsizes(Times)
rank = dimsizes(dimT)
if (rank.ne.2) then
print("===> WRF_contributed.ncl: WRF_Times2Udunits_c expects 2D array: rank="+rank)
exit
end if
if (.not.(typeof(opt).eq."integer" .or. typeof(opt).eq."string")) then
print("===> WRF_contributed.ncl: opt must be integer or string: type="+typeof(opt))
exit
end if
year = stringtointeger((/Times(:, 0:3) /))
month = stringtointeger((/Times(:, 5:6) /))
day = stringtointeger((/Times(:, 8:9) /))
hour = stringtointeger((/Times(:,11:12)/))
minute = stringtointeger((/Times(:,14:15)/))
sec = stringtointeger((/Times(:,17:18)/))
if (typeof(opt).eq."integer") then
if (opt.eq.0) then
units = "hours since "+year(0)+"-" \
+sprinti("%0.2i",month(0)) +"-" \
+sprinti("%0.2i",day(0)) +" " \
+sprinti("%0.2i",hour(0)) +":" \
+sprinti("%0.2i",minute(0))+":" \
+sprinti("%0.2i",sec(0))
else
units = "hours since 1901-01-01 00:00:00"
end if
else
units = opt ; opt is udunits compatible string
end if
Time = ut_inv_calendar(year,month,day,hour,minute,sec, units, 0)
Time!0 = "Time"
Time@long_name = "Time"
Time@description= "Time"
Time@units = units
Time&Time = Time ; make coordinate variable
return (Time)
end
;*************************************************************************
; D. Shea
; interface to WRF_Times2Udunits_c
undef ("WRF_Times_to_udunits")
function WRF_Times_to_udunits(Times:character, opt)
begin
return( WRF_Times2Udunits_c(Times, 0) ) ; old name
end
;*************************************************************************
; D. Shea
; convert WRF character variable "Times" to
; a coordinate variable of type double
; time(double) = yyyymmddhhmnss
; 2001-06-11_12:00:00 ==> 20010611120000
;
; opt: currently not used [dummy]
;
undef ("WRF_Times2double_c")
function WRF_Times2double_c(Times:character, opt)
local dimT, rank, N, time, i, tmp_c
begin
dimT = dimsizes(Times)
rank = dimsizes(dimT)
if (rank.ne.2) then
print("===> WRF_contributed.ncl: WRF_Times2Udunits_c expects 2D array: rank="+rank)
exit
end if
N = dimT(0)
Time = new( N ,"double") ; preset to "double"
delete(Time@_FillValue) ; coord variables should not have a _FillValue
Time = stringtointeger((/Times(:,0:3)/)) *1d10 + \ ; yyyy
stringtointeger((/Times(:,5:6)/)) *1d8 + \ ; mm
stringtointeger((/Times(:,8:9)/)) *1d6 + \ ; dd
stringtointeger((/Times(:,11:12)/))*1d4 + \ ; hh
stringtointeger((/Times(:,14:15)/))*1d2 + \ ; mn
stringtointeger((/Times(:,17:18)/))*1d0 ; ss
Time!0 = "Time"
Time@long_name = "Time"
Time@description= "Time"
Time@units = "yyyymmddhhmnss"
Time&Time = Time ; make coordinate variable
return (Time)
end
;*************************************************************************
; D. Shea
; interface to WRF_Times2double_c
; more explicit function name
undef ("WRF_Times_to_ymdhms")
function WRF_Times_to_ymdhms(Times:character, opt)
begin
return( WRF_Times2double_c(Times, 0) ) ; old name
end
;*************************************************************************
; D. Shea
; convert WRF character variable "Times" to
; a coordinate variable of type integer
; time(integer)= yyyymmddhh [->ymdh]
; 2001-06-11_12:00:00 ==> 2001061112
;
; Note: mminute and second are not part of the returned time
;
; opt: currently not used [dummy]
;
undef ("WRF_Times_to_ymdh")
function WRF_Times_to_ymdh(Times:character, opt)
local dimT, rank, N, time, i, tmp_c
begin
dimT = dimsizes(Times)
rank = dimsizes(dimT)
if (rank.ne.2) then
print("===> WRF_contributed.ncl: WRF_Times2yyyymmddhh expects 2D array: rank="+rank)
exit
end if
N = dimT(0)
Time = new( N ,"integer")
delete(Time@_FillValue) ; coord variables should not have a _FillValue
Time = stringtointeger((/Times(:,0:3)/)) *1000000 + \ ; yyyy
stringtointeger((/Times(:,5:6)/)) *10000 + \ ; mm
stringtointeger((/Times(:,8:9)/)) *100 + \ ; dd
stringtointeger((/Times(:,11:12)/)) ; hh
Time!0 = "Time"
Time@long_name = "Time"
Time@description= "Time"
Time@units = "yyyymmddhh"
Time&Time = Time ; make coordinate variable
return (Time)
end
;*************************************************************************
; D. Shea
; This is a driver that selects the appropriate
; mapping function based upon the file attribute: MAP_PROJ
; MAP_PROJ=1 [Lambert Conformal]; =2 [Stereographic]; =3 [Mercator]
;
; opt: currently not used [potentail use: time counter for XLAT/XLONG]
;
; Sample usage:
; ncdf = addfile("...", r")
; res = True
; WRF_map_c (ncdf, res, 0)
; res = ...
;
undef("WRF_map_c")
procedure WRF_map_c (f:file, res:logical, opt)
local rank, dimll, nlat, mlon, lat2d, lon2d
begin
if (isatt(f,"MAP_PROJ")) then
if (f@MAP_PROJ.eq.1) then
res@mpProjection = "LambertConformal"
end if
if (f@MAP_PROJ.eq.2) then
res@mpProjection = "Stereographic"
end if
if (f@MAP_PROJ.eq.3) then
res@mpProjection = "Mercator"
end if
else
print ("WRF_mapProj: no MAP_PROJ attribute")
end if
rank = dimsizes(filevardimsizes(f,"XLAT")) ; # of dimensions
if (rank.eq.3) then
lat2d = f->XLAT(0,:,:) ; opt could bt "nt" f->XLAT(opt,:,:)
lon2d = f->XLONG(0,:,:)
else
if (rank.eq.2) then
lat2d = f->XLAT
lon2d = f->XLONG
else
print ("WRF_resLamCon_c: unexpected lat/lon rank: rank="+rank)
exit
end if
end if
lat2d@units = "degrees_north" ; not needed
lon2d@units = "degrees_east"
dimll = dimsizes(lat2d)
nlat = dimll(0)
mlon = dimll(1)
res@mpLimitMode = "Corners"
res@mpLeftCornerLatF = lat2d(0,0)
res@mpLeftCornerLonF = lon2d(0,0)
res@mpRightCornerLatF = lat2d(nlat-1,mlon-1)
res@mpRightCornerLonF = lon2d(nlat-1,mlon-1)
res@mpCenterLonF = f@CEN_LON
res@mpCenterLatF = f@CEN_LAT ; default
if (res@mpProjection.eq."Mercator") then
res@mpCenterLatF = 0.0 ; Cindy Bruyere MMM/WRF 24 Mar 2006
end if
if (res@mpProjection.eq."LambertConformal") then
res@mpLambertParallel1F = f@TRUELAT1
res@mpLambertParallel2F = f@TRUELAT2
if (isatt(f, "STAND_LON") ) then
res@mpLambertMeridianF = f@STAND_LON ; use if present
; CB MMM/WRF 4 Aug 2006
else
if (isatt(f, "CEN_LON") ) then
res@mpLambertMeridianF = f@CEN_LON
else
print("WRF_map_c: STAND_LON and CEN_LON missing")
end if
end if
end if
res@mpFillOn = False ; turn off map fill
res@mpOutlineDrawOrder = "PostDraw" ; draw continental outline last
res@mpOutlineBoundarySets = "GeophysicalAndUSStates" ; state boundaries
res@mpPerimDrawOrder = "PostDraw" ; force map perim
; commented 5/17/2007
;;res@tfDoNDCOverlay = True ; True for 'native' grid
; some WRF are not native
res@gsnAddCyclic = False ; data are not cyclic
end
;*************************************************************************
; D. Shea
; interface for backward compatibility
undef("WRF_resLamCon_c")
procedure WRF_resLamCon_c (f:file, res:logical, opt)
begin
WRF_map_c (f, res, opt)
end
;*************************************************************************
; D. Shea
; interface for newly named procedure
undef("wrf_mapres_c")
procedure wrf_mapres_c(f:file, res:logical, opt)
begin
WRF_map_c (f, res, opt)
end
;*************************************************************************
; D. Shea
; single interface to convert WRF character variable "Times"
; to user specified numeric values
;
; M. Haley
; At some point we decided to rename this from WRF_Times to wrf_times_c
; Also added error check for opt.
;
undef ("wrf_times_c")
function wrf_times_c(Times:character, opt:integer)
begin
if (opt.ge.0 .and. opt.le.1) then
return(WRF_Times2Udunits_c(Times, opt) )
end if
if (opt.eq.2) then
return(WRF_Times2double_c(Times, opt) )
end if
if (opt.eq.3) then
return(WRF_Times_to_ymdh(Times, opt) )
end if
end

158
ncl_reference/calc_uh.f90
Unescape View File

@ -1,158 +0,0 @@
! For NCL graphics:
! WRAPIT -m64 calc_uh90.stub calc_uh.f90
! This should create a shared library named "calc_uh90.so".
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE CALC_UH ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
! Calculates updraft helicity (UH) to detect rotating updrafts.
! Formula follows Kain et al, 2008, Wea. and Forecasting, 931-952,
! but this version has controls for the limits of integration
! uhminhgt to uhmxhgt, in m AGL. Kain et al used 2000 to 5000 m.
! Units of UH are m^2/s^2.
!
! Note here that us and vs are at ARPS scalar points.
! w is at w-point (scalar pt in horiz, staggered vertical)
!
! Keith Brewster, CAPS/Univ. of Oklahoma
! March, 2010
!
! uh = wrf_updraft_helicity(zp,us,vs,w,
SUBROUTINE dcalcuh(nx,ny,nz,nzp1,zp,mapfct,dx,dy,uhmnhgt,uhmxhgt, &
us,vs,w,uh,tem1,tem2)
IMPLICIT NONE
INTEGER, INTENT(IN) :: nx,ny,nz,nzp1
DOUBLE PRECISION, INTENT(IN) :: zp(nx,ny,nzp1)
DOUBLE PRECISION, INTENT(IN) :: mapfct(nx,ny)
DOUBLE PRECISION, INTENT(IN) :: dx,dy
DOUBLE PRECISION, INTENT(IN) :: uhmnhgt,uhmxhgt
DOUBLE PRECISION, INTENT(IN) :: us(nx,ny,nz)
DOUBLE PRECISION, INTENT(IN) :: vs(nx,ny,nz)
DOUBLE PRECISION, INTENT(IN) :: w(nx,ny,nzp1)
DOUBLE PRECISION, INTENT(OUT) :: uh(nx,ny)
DOUBLE PRECISION, INTENT(OUT) :: tem1(nx,ny,nz)
DOUBLE PRECISION, INTENT(OUT) :: tem2(nx,ny,nz)
!
! Misc local variables
!
INTEGER :: i,j,k,kbot,ktop
DOUBLE PRECISION :: twodx,twody,wgtlw,sum,wmean,wsum,wavg
DOUBLE PRECISION :: helbot,heltop,wbot,wtop
DOUBLE PRECISION :: zbot,ztop
!
! Initialize arrays
!
uh=0.0
tem1=0.0
!
! Calculate vertical component of helicity at scalar points
! us: u at scalar points
! vs: v at scalar points
!
twodx=2.0*dx
twody=2.0*dy
DO k=2,nz-2
DO j=2,ny-1
DO i=2,nx-1
wavg=0.5*(w(i,j,k)+w(i,j,k+1))
tem1(i,j,k)=wavg * &
((vs(i+1,j,k)-vs(i-1,j,k))/(twodx*mapfct(i,j)) - &
(us(i,j+1,k)-us(i,j-1,k))/(twody*mapfct(i,j)))
tem2(i,j,k)=0.5*(zp(i,j,k)+zp(i,j,k+1))
END DO
END DO
END DO
!
! Integrate over depth uhminhgt to uhmxhgt AGL
!
! WRITE(6,'(a,f12.1,a,f12.1,a)') &
! 'Calculating UH from ',uhmnhgt,' to ',uhmxhgt,' m AGL'
DO j=2,ny-2
DO i=2,nx-2
zbot=zp(i,j,2)+uhmnhgt
ztop=zp(i,j,2)+uhmxhgt
!
! Find wbar, weighted-mean vertical velocity in column
! Find w at uhmnhgt AGL (bottom)
!
DO k=2,nz-3
IF(zp(i,j,k) > zbot) EXIT
END DO
kbot=k
wgtlw=(zp(i,j,kbot)-zbot)/(zp(i,j,kbot)-zp(i,j,kbot-1))
wbot=(wgtlw*w(i,j,kbot-1))+((1.-wgtlw)*w(i,j,kbot))
!
! Find w at uhmxhgt AGL (top)
!
DO k=2,nz-3
IF(zp(i,j,k) > ztop) EXIT
END DO
ktop=k
wgtlw=(zp(i,j,ktop)-ztop)/(zp(i,j,ktop)-zp(i,j,ktop-1))
wtop=(wgtlw*w(i,j,ktop-1))+((1.-wgtlw)*w(i,j,ktop))
!
! First part, uhmnhgt to kbot
!
wsum=0.5*(w(i,j,kbot)+wbot)*(zp(i,j,kbot)-zbot)
!
! Integrate up through column
!
DO k=(kbot+1),(ktop-1)
wsum=wsum+0.5*(w(i,j,k)+w(i,j,k-1))*(zp(i,j,k)-zp(i,j,k-1))
END DO
!
! Last part, ktop-1 to uhmxhgt
!
wsum=wsum+0.5*(wtop+w(i,j,ktop-1))*(ztop-zp(i,j,ktop-1))
wmean=wsum/(uhmxhgt-uhmnhgt)
IF(wmean > 0.) THEN ! column updraft, not downdraft
!
! Find helicity at uhmnhgt AGL (bottom)
!
DO k=2,nz-3
IF(tem2(i,j,k) > zbot) EXIT
END DO
kbot=k
wgtlw=(tem2(i,j,kbot)-zbot)/(tem2(i,j,kbot)-tem2(i,j,kbot-1))
helbot=(wgtlw*tem1(i,j,kbot-1))+((1.-wgtlw)*tem1(i,j,kbot))
!
! Find helicity at uhmxhgt AGL (top)
!
DO k=2,nz-3
IF(tem2(i,j,k) > ztop) EXIT
END DO
ktop=k
wgtlw=(tem2(i,j,ktop)-ztop)/(tem2(i,j,ktop)-tem2(i,j,ktop-1))
heltop=(wgtlw*tem1(i,j,ktop-1))+((1.-wgtlw)*tem1(i,j,ktop))
!
! First part, uhmnhgt to kbot
!
sum=0.5*(tem1(i,j,kbot)+helbot)*(tem2(i,j,kbot)-zbot)
!
! Integrate up through column
!
DO k=(kbot+1),(ktop-1)
sum=sum+0.5*(tem1(i,j,k)+tem1(i,j,k-1))*(tem2(i,j,k)-tem2(i,j,k-1))
END DO
!
! Last part, ktop-1 to uhmxhgt
!
uh(i,j)=sum+0.5*(heltop+tem1(i,j,ktop-1))*(ztop-tem2(i,j,ktop-1))
END IF
END DO
END DO
uh = uh * 1000. ! Scale according to Kain et al. (2008)
RETURN
END SUBROUTINE dcalcuh

60
ncl_reference/cloud_fracf.f
Unescape View File

@ -1,60 +0,0 @@
C NCLFORTSTART
subroutine cloud_frac(pres,rh,lowc,midc,highc,nz,ns,ew)
implicit none
integer nz,ns,ew
real pres(ew,ns,nz),rh(ew,ns,nz)
real lowc(ew,ns),midc(ew,ns),highc(ew,ns)
C NCLEND
integer i,j,k
integer kchi,kcmi,kclo
DO j = 1,ns
DO i = 1,ew
DO k = 1,nz-1
c if((pres(i,j,k) .ge. 45000. ) .and.
c & (pres(i,j,k) .lt. 80000.)) then
c kchi = k
c else if((pres(i,j,k) .ge. 80000.) .and.
c & (pres(i,j,k) .lt. 97000.)) then
c kcmi = k
c else if (pres(i,j,k) .ge. 97000.) then
c kclo = k
c end if
IF ( pres(i,j,k) .gt. 97000. ) kclo=k
IF ( pres(i,j,k) .gt. 80000. ) kcmi=k
IF ( pres(i,j,k) .gt. 45000. ) kchi=k
end do
DO k = 1,nz-1
IF ( k .ge. kclo .AND. k .lt. kcmi ) then
lowc(i,j) = AMAX1(rh(i,j,k),lowc(i,j))
else IF ( k .ge. kcmi .AND. k .lt. kchi ) then !! mid cloud
midc(i,j) = AMAX1(rh(i,j,k),midc(i,j))
else if ( k .ge. kchi ) then !! high cloud
highc(i,j) = AMAX1(rh(i,j,k),highc(i,j))
end if
END DO
lowc(i,j) = 4.0 * lowc(i,j)/100.-3.0
midc(i,j) = 4.0 * midc(i,j)/100.-3.0
highc(i,j) = 2.5 * highc(i,j)/100.-1.5
lowc(i,j) = amin1(lowc(i,j),1.0)
lowc(i,j) = amax1(lowc(i,j),0.0)
midc(i,j) = amin1(midc(i,j),1.0)
midc(i,j) = amax1(midc(i,j),0.0)
highc(i,j) = amin1(highc(i,j),1.0)
highc(i,j) = amax1(highc(i,j),0.0)
END DO
END DO
return
end

72
ncl_reference/eqthecalc.f
Unescape View File

@ -1,72 +0,0 @@
SUBROUTINE DEQTHECALC(QVP,TMK,PRS,ETH,MIY,MJX,MKZH)
DOUBLE PRECISION EPS
DOUBLE PRECISION RGAS
DOUBLE PRECISION RGASMD
DOUBLE PRECISION CP
DOUBLE PRECISION CPMD
DOUBLE PRECISION GAMMA
DOUBLE PRECISION GAMMAMD
DOUBLE PRECISION TLCLC1
DOUBLE PRECISION TLCLC2
DOUBLE PRECISION TLCLC3
DOUBLE PRECISION TLCLC4
DOUBLE PRECISION THTECON1
DOUBLE PRECISION THTECON2
DOUBLE PRECISION THTECON3
DOUBLE PRECISION Q
DOUBLE PRECISION T
DOUBLE PRECISION P
DOUBLE PRECISION E
DOUBLE PRECISION TLCL
c
c Input variables
c Qvapor [g/kg]
DOUBLE PRECISION QVP(MIY,MJX,MKZH)
c Temperature [K]
DOUBLE PRECISION TMK(MIY,MJX,MKZH)
c full pressure (=P+PB) [hPa]
DOUBLE PRECISION PRS(MIY,MJX,MKZH)
c
c Output variable
c equivalent potential temperature [K]
DOUBLE PRECISION ETH(MIY,MJX,MKZH)
c
c parameters
PARAMETER (EPS=0.622D0)
c J/K/kg
RGAS = 287.04D0
c rgas_moist=rgas*(1.+rgasmd*qvp)
RGASMD = .608D0
c J/K/kg Note: not using Bolton's value of 1005.7
CP = 1004.D0
c cp_moist=cp*(1.+cpmd*qvp)
CPMD = .887D0
GAMMA = RGAS/CP
c gamma_moist=gamma*(1.+gammamd*qvp)
GAMMAMD = RGASMD - CPMD
TLCLC1 = 2840.D0
TLCLC2 = 3.5D0
TLCLC3 = 4.805D0
TLCLC4 = 55.D0
c K
THTECON1 = 3376.D0
THTECON2 = 2.54D0
THTECON3 = .81D0
c
DO 1000 K = 1,MKZH
DO 1000 J = 1,MJX
DO 1000 I = 1,MIY
Q = MAX(QVP(I,J,K),1.D-15)
T = TMK(I,J,K)
P = PRS(I,J,K)/100.
E = Q*P/ (EPS+Q)
TLCL = TLCLC1/ (LOG(T**TLCLC2/E)-TLCLC3) + TLCLC4
ETH(I,J,K) = T* (1000.D0/P)**
+ (GAMMA* (1.D0+GAMMAMD*Q))*
+ EXP((THTECON1/TLCL-THTECON2)*Q*
+ (1.D0+THTECON3*Q))
1000 CONTINUE
RETURN
END

82
ncl_reference/etaconv.py
Unescape View File

@ -1,82 +0,0 @@
import numpy as n
from wrf.var.extension import computeeta
from wrf.var.constants import Constants
from wrf.var.decorators import convert_units
from wrf.var.util import extract_vars
#__all__ = ["convert_eta"]
__all__ = []
# A useful utility, but should probably just use geopotential height when
# plotting for AGL levels
# Eta definition (nu):
# nu = (P - Ptop) / (Psfc - Ptop)
# def convert_eta(wrfnc, p_or_z="ht", timeidx=0):
# if (p_or_z.lower() == "height" or p_or_z.lower() == "ht"
# or p_or_z.lower() == "h"):
# return_z = True
# elif (p_or_z.lower() == "p" or p_or_z.lower() == "pres"
# or p_or_z.lower() == "pressure"):
# return_z = False
#
# R = Constants.R
# G = Constants.G
# CP = Constants.CP
#
# # Keeping the slice notation to show the dimensions
# # Note: Not sure if T00 should be used (290) or the usual hard-coded 300 for base
# # theta
# height_data = wrfnc.variables["HGT"][timeidx,:,:]
# znu_data = wrfnc.variables["ZNU"][timeidx,:]
# #t00_data = wrfnc.variables["T00"][timeidx]
# psfc_data = wrfnc.variables["PSFC"][timeidx,:,:]
# ptop_data = wrfnc.variables["P_TOP"][timeidx]
# pth_data = wrfnc.variables["T"][timeidx,:,:,:] # Pert potential temp
#
# pcalc_data = n.zeros(pth_data.shape, dtype=n.float32)
# mean_t_data = n.zeros(pth_data.shape, dtype=n.float32)
# temp_data = n.zeros(pth_data.shape, dtype=n.float32)
# z_data = n.zeros(pth_data.shape, dtype=n.float32)
#
# #theta_data = pth_data + t00_data
# theta_data = pth_data + Constants.T_BASE
#
# for k in xrange(znu_data.shape[0]):
# pcalc_data[k,:,:] = znu_data[k]*(psfc_data[:,:] - (ptop_data)) + (ptop_data)
#
# # Potential temperature:
# # theta = T * (Po/P)^(R/CP)
# #
# # Hypsometric equation:
# # h = z2-z1 = R*Tbar/G * ln(p1/p2)
# # where z1, p1 are the surface
# if return_z:
# for k in xrange(znu_data.shape[0]):
# temp_data[k,:,:] = (theta_data[k,:,:]) / ((100000.0 / (pcalc_data[k,:,:]))**(R/CP))
# mean_t_data[k,:,:] = n.mean(temp_data[0:k+1,:,:], axis=0)
# z_data[k,:,:] = ((R*mean_t_data[k,:,:]/G) * n.log(psfc_data[:,:]/pcalc_data[k,:,:]))
#
# return z_data
# else:
# return pcalc_data * .01
# def convert_eta(wrfnc, units="m", msl=False, timeidx=0):
# check_units(units, "height")
# hgt = wrfnc.variables["HGT"][timeidx,:,:]
# znu = wrfnc.variables["ZNU"][timeidx,:]
# psfc = wrfnc.variables["PSFC"][timeidx,:,:]
# ptop = wrfnc.variables["P_TOP"][timeidx]
# t = wrfnc.variables["T"][timeidx,:,:,:]
#
# full_theta = t + Constants.T_BASE
#
# z = computeeta(full_theta, znu, psfc, ptop)
#
# if not msl:
# return convert_units(z, "height", "m", units)
# else:
# return convert_units(z + hgt, "height", "m", units)

183
ncl_reference/jinja_context.py
Unescape View File

@ -1,183 +0,0 @@
'''
Created on Jan 16, 2014
@author: sean
'''
from __future__ import absolute_import, division, print_function
from functools import partial
import json
import logging
import os
import sys
import jinja2
from .conda_interface import PY3
from .environ import get_dict as get_environ
from .metadata import select_lines, ns_cfg
from .source import WORK_DIR
log = logging.getLogger(__file__)
class UndefinedNeverFail(jinja2.Undefined):
"""
A class for Undefined jinja variables.
This is even less strict than the default jinja2.Undefined class,
because it permits things like {{ MY_UNDEFINED_VAR[:2] }} and
{{ MY_UNDEFINED_VAR|int }}. This can mask lots of errors in jinja templates, so it
should only be used for a first-pass parse, when you plan on running a 'strict'
second pass later.
"""
all_undefined_names = []
def __init__(self, hint=None, obj=jinja2.runtime.missing, name=None,
exc=jinja2.exceptions.UndefinedError):
UndefinedNeverFail.all_undefined_names.append(name)
jinja2.Undefined.__init__(self, hint, obj, name, exc)
__add__ = __radd__ = __mul__ = __rmul__ = __div__ = __rdiv__ = \
__truediv__ = __rtruediv__ = __floordiv__ = __rfloordiv__ = \
__mod__ = __rmod__ = __pos__ = __neg__ = __call__ = \
__getitem__ = __lt__ = __le__ = __gt__ = __ge__ = \
__complex__ = __pow__ = __rpow__ = \
lambda self, *args, **kwargs: UndefinedNeverFail(hint=self._undefined_hint,
obj=self._undefined_obj,
name=self._undefined_name,
exc=self._undefined_exception)
__str__ = __repr__ = \
lambda *args, **kwargs: u''
__int__ = lambda _: 0
__float__ = lambda _: 0.0
def __getattr__(self, k):
try:
return object.__getattr__(self, k)
except AttributeError:
return UndefinedNeverFail(hint=self._undefined_hint,
obj=self._undefined_obj,
name=self._undefined_name + '.' + k,
exc=self._undefined_exception)
class FilteredLoader(jinja2.BaseLoader):
"""
A pass-through for the given loader, except that the loaded source is
filtered according to any metadata selectors in the source text.
"""
def __init__(self, unfiltered_loader):
self._unfiltered_loader = unfiltered_loader
self.list_templates = unfiltered_loader.list_templates
def get_source(self, environment, template):
contents, filename, uptodate = self._unfiltered_loader.get_source(environment,
template)
return select_lines(contents, ns_cfg()), filename, uptodate
def load_setup_py_data(setup_file='setup.py', from_recipe_dir=False, recipe_dir=None,
unload_modules=None, fail_on_error=False):
_setuptools_data = {}
def setup(**kw):
_setuptools_data.update(kw)
import setuptools
import distutils.core
try:
import numpy.distutils.core
except ImportError:
do_numpy = False
else:
do_numpy = True
cd_to_work = False
if from_recipe_dir and recipe_dir:
setup_file = os.path.abspath(os.path.join(recipe_dir, setup_file))
elif os.path.exists(WORK_DIR):
cd_to_work = True
cwd = os.getcwd()
os.chdir(WORK_DIR)
if not os.path.isabs(setup_file):
setup_file = os.path.join(WORK_DIR, setup_file)
# this is very important - or else if versioneer or otherwise is in the start folder,
# things will pick up the wrong versioneer/whatever!
sys.path.insert(0, WORK_DIR)
else:
log.debug("Did not find setup.py file in manually specified location, and source "
"not downloaded yet.")
return {}
# Patch setuptools, distutils
setuptools_setup = setuptools.setup
distutils_setup = distutils.core.setup
setuptools.setup = distutils.core.setup = setup
if do_numpy:
numpy_setup = numpy.distutils.core.setup
numpy.distutils.core.setup = setup
ns = {
'__name__': '__main__',
'__doc__': None,
'__file__': setup_file,
}
try:
code = compile(open(setup_file).read(), setup_file, 'exec', dont_inherit=1)
exec(code, ns, ns)
distutils.core.setup = distutils_setup
setuptools.setup = setuptools_setup
if do_numpy:
numpy.distutils.core.setup = numpy_setup
# this happens if setup.py is used in load_setup_py_data, but source is not yet downloaded
except:
raise
finally:
if cd_to_work:
os.chdir(cwd)
del sys.path[0]
return _setuptools_data
def load_setuptools(setup_file='setup.py', from_recipe_dir=False, recipe_dir=None,
unload_modules=None, fail_on_error=False):
log.warn("Deprecation notice: the load_setuptools function has been renamed to "
"load_setup_py_data. load_setuptools will be removed in a future release.")
return load_setup_py_data(setup_file=setup_file, from_recipe_dir=from_recipe_dir,
recipe_dir=recipe_dir, unload_modules=unload_modules,
fail_on_error=fail_on_error)
def load_npm():
# json module expects bytes in Python 2 and str in Python 3.
mode_dict = {'mode': 'r', 'encoding': 'utf-8'} if PY3 else {'mode': 'rb'}
with open('package.json', **mode_dict) as pkg:
return json.load(pkg)
def context_processor(initial_metadata, recipe_dir):
"""
Return a dictionary to use as context for jinja templates.
initial_metadata: Augment the context with values from this MetaData object.
Used to bootstrap metadata contents via multiple parsing passes.
"""
ctx = get_environ(m=initial_metadata)
environ = dict(os.environ)
environ.update(get_environ(m=initial_metadata))
ctx.update(
load_setup_py_data=partial(load_setup_py_data, recipe_dir=recipe_dir),
# maintain old alias for backwards compatibility:
load_setuptools=partial(load_setuptools, recipe_dir=recipe_dir),
load_npm=load_npm,
environ=environ)
return ctx

170
ncl_reference/module_model_constants.f
Unescape View File

@ -1,170 +0,0 @@
!WRF:MODEL_LAYER:CONSTANTS
!
MODULE module_model_constants
! 2. Following are constants for use in defining real number bounds.
! A really small number.
REAL , PARAMETER :: epsilon = 1.E-15
! 4. Following is information related to the physical constants.
! These are the physical constants used within the model.
! JM NOTE -- can we name this grav instead?
REAL , PARAMETER :: g = 9.81 ! acceleration due to gravity (m {s}^-2)
#if ( NMM_CORE == 1 )
REAL , PARAMETER :: r_d = 287.04
REAL , PARAMETER :: cp = 1004.6
#else
REAL , PARAMETER :: r_d = 287.
REAL , PARAMETER :: cp = 7.*r_d/2.
#endif
REAL , PARAMETER :: r_v = 461.6
REAL , PARAMETER :: cv = cp-r_d
REAL , PARAMETER :: cpv = 4.*r_v
REAL , PARAMETER :: cvv = cpv-r_v
REAL , PARAMETER :: cvpm = -cv/cp
REAL , PARAMETER :: cliq = 4190.
REAL , PARAMETER :: cice = 2106.
REAL , PARAMETER :: psat = 610.78
REAL , PARAMETER :: rcv = r_d/cv
REAL , PARAMETER :: rcp = r_d/cp
REAL , PARAMETER :: rovg = r_d/g
REAL , PARAMETER :: c2 = cp * rcv
real , parameter :: mwdry = 28.966 ! molecular weight of dry air (g/mole)
REAL , PARAMETER :: p1000mb = 100000.
REAL , PARAMETER :: t0 = 300.
REAL , PARAMETER :: p0 = p1000mb
REAL , PARAMETER :: cpovcv = cp/(cp-r_d)
REAL , PARAMETER :: cvovcp = 1./cpovcv
REAL , PARAMETER :: rvovrd = r_v/r_d
REAL , PARAMETER :: reradius = 1./6370.0e03
REAL , PARAMETER :: asselin = .025
! REAL , PARAMETER :: asselin = .0
REAL , PARAMETER :: cb = 25.
REAL , PARAMETER :: XLV0 = 3.15E6
REAL , PARAMETER :: XLV1 = 2370.
REAL , PARAMETER :: XLS0 = 2.905E6
REAL , PARAMETER :: XLS1 = 259.532
REAL , PARAMETER :: XLS = 2.85E6
REAL , PARAMETER :: XLV = 2.5E6
REAL , PARAMETER :: XLF = 3.50E5
REAL , PARAMETER :: rhowater = 1000.
REAL , PARAMETER :: rhosnow = 100.
REAL , PARAMETER :: rhoair0 = 1.28
!
! Now namelist-specified parameter: ccn_conc - RAS
! REAL , PARAMETER :: n_ccn0 = 1.0E8
!
REAL , PARAMETER :: piconst = 3.1415926535897932384626433
REAL , PARAMETER :: DEGRAD = piconst/180.
REAL , PARAMETER :: DPD = 360./365.
REAL , PARAMETER :: SVP1=0.6112
REAL , PARAMETER :: SVP2=17.67
REAL , PARAMETER :: SVP3=29.65
REAL , PARAMETER :: SVPT0=273.15
REAL , PARAMETER :: EP_1=R_v/R_d-1.
REAL , PARAMETER :: EP_2=R_d/R_v
REAL , PARAMETER :: KARMAN=0.4
REAL , PARAMETER :: EOMEG=7.2921E-5
REAL , PARAMETER :: STBOLT=5.67051E-8
REAL , PARAMETER :: prandtl = 1./3.0
! constants for w-damping option
REAL , PARAMETER :: w_alpha = 0.3 ! strength m/s/s
REAL , PARAMETER :: w_beta = 1.0 ! activation cfl number
REAL , PARAMETER :: pq0=379.90516
REAL , PARAMETER :: epsq2=0.2
REAL , PARAMETER :: a2=17.2693882
REAL , PARAMETER :: a3=273.16
REAL , PARAMETER :: a4=35.86
REAL , PARAMETER :: epsq=1.e-12
REAL , PARAMETER :: p608=rvovrd-1.
!#if ( NMM_CORE == 1 )
REAL , PARAMETER :: climit=1.e-20
REAL , PARAMETER :: cm1=2937.4
REAL , PARAMETER :: cm2=4.9283
REAL , PARAMETER :: cm3=23.5518
! REAL , PARAMETER :: defc=8.0
! REAL , PARAMETER :: defm=32.0
REAL , PARAMETER :: defc=0.0
REAL , PARAMETER :: defm=99999.0
REAL , PARAMETER :: epsfc=1./1.05
REAL , PARAMETER :: epswet=0.0
REAL , PARAMETER :: fcdif=1./3.
#if ( HWRF == 1 )
REAL , PARAMETER :: fcm=0.0
#else
REAL , PARAMETER :: fcm=0.00003
#endif
REAL , PARAMETER :: gma=-r_d*(1.-rcp)*0.5
REAL , PARAMETER :: p400=40000.0
REAL , PARAMETER :: phitp=15000.0
REAL , PARAMETER :: pi2=2.*3.1415926, pi1=3.1415926
REAL , PARAMETER :: plbtm=105000.0
REAL , PARAMETER :: plomd=64200.0
REAL , PARAMETER :: pmdhi=35000.0
REAL , PARAMETER :: q2ini=0.50
REAL , PARAMETER :: rfcp=0.25/cp
REAL , PARAMETER :: rhcrit_land=0.75
REAL , PARAMETER :: rhcrit_sea=0.80
REAL , PARAMETER :: rlag=14.8125
REAL , PARAMETER :: rlx=0.90
REAL , PARAMETER :: scq2=50.0
REAL , PARAMETER :: slopht=0.001
REAL , PARAMETER :: tlc=2.*0.703972477
REAL , PARAMETER :: wa=0.15
REAL , PARAMETER :: wght=0.35
REAL , PARAMETER :: wpc=0.075
REAL , PARAMETER :: z0land=0.10
#if ( HWRF == 1 )
REAL , PARAMETER :: z0max=0.01
#else
REAL , PARAMETER :: z0max=0.008
#endif
REAL , PARAMETER :: z0sea=0.001
!#endif
! Earth
! The value for P2SI *must* be set to 1.0 for Earth
! Although, now we may not need this declaration here (see above)
!REAL , PARAMETER :: P2SI = 1.0
! Orbital constants:
INTEGER , PARAMETER :: PLANET_YEAR = 365
REAL , PARAMETER :: OBLIQUITY = 23.5
REAL , PARAMETER :: ECCENTRICITY = 0.014
REAL , PARAMETER :: SEMIMAJORAXIS = 1.0 ! In AU
! Don't know the following values, so we'll fake them for now
REAL , PARAMETER :: zero_date = 0.0 ! Time of perihelion passage
! Fraction into the year (from perhelion) of the
! occurrence of the Northern Spring Equinox
REAL , PARAMETER :: EQUINOX_FRACTION= 0.0
! 2012103
#if (EM_CORE == 1)
! for calls to set_tiles
INTEGER, PARAMETER :: ZONE_SOLVE_EM = 1
INTEGER, PARAMETER :: ZONE_SFS = 2
#endif
CONTAINS
SUBROUTINE init_module_model_constants
END SUBROUTINE init_module_model_constants
END MODULE module_model_constants

77
ncl_reference/product_table.txt
Unescape View File

@ -1,77 +0,0 @@
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| Variable Name | Description | Units | Additional Keyword Arguments |
+====================+===============================================================+=====================+===============================================================================================+
| avo | Absolute Vorticity | 10-5 s-1 | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| eth/theta_e | Equivalent Potential Temperature | K | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| cape_2d | 2D cape (mcape/mcin/lcl/lfc) | J/kg / J/kg / m / m | missing: Fill value for output only (float) |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| cape_3d | 3D cape and cin | J/kg | missing: Fill value for output only (float) |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| ctt | Cloud Top Temperature | C | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| cloudfrac | Cloud Fraction | % | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| dbz | Reflectivity | dBz | do_variant: Set to True to enable variant calculation. Default is False. |
| | | | do_liqskin : Set to True to enable liquid skin calculation. Default is False. |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| mdbz | Maximum Reflectivity | dBz | do_variant: Set to True to enable variant calculation. Default is False. |
| | | | do_liqskin: Set to True to enable liquid skin calculation. Default is False. |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| geopt/geopotential | Full Model Geopotential | m2 s-2 | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| helicity | Storm Relative Helicity | m2 s-2 | top: The top level for the calculation in meters (float). Default is 3000.0. |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| lat | Latitude | decimal degrees | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| lon | Longitude | decimal degrees | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| omg/omega | Omega | Pa/s | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| p/pres | Full Model Pressure | Pa | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| pressure | Full Model Pressure | hPa | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| pvo | Potential Vorticity | PVU | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| pw | Precipitable Water | kg m-2 | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| rh2 | 2m Relative Humidity | % | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| slp | Sea Level Pressure | hPa | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| ter | Model Terrain Height | m | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| td2 | 2m Dew Point Temperature | C | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| td | Dew Point Temperature | C | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| tc | Temperature | C | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| th/theta | Potential Temperature | K | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| tk | Temperature | K | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| times | Times in the File or Sequence | | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| tv | Virtual Temperature | K | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| twb | Wet Bulb Temperature | K | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| updraft_helicity | Updraft Helicity | m2 s-2 | bottom: The bottom level for the calculation in meters (float). Default is 2000.0. |
| | | | top: The top level for the calculation in meters (float). Default is 5000.0. |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| ua | U-component of Wind on Mass Points | m/s | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| va | V-component of Wind on Mass Points | m/s | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| wa | W-component of Wind on Mass Points | m/s | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| uvmet10 | 10 m U and V Components of Wind Rotated to Earth Coordinates | m/s | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| uvmet | U and V Components of Wind Rotated to Earth Coordinates | m/s | |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+
| z/height | Full Model Height | m | msl: Set to False to return AGL values. Otherwise, MSL. Default is True. |
+--------------------+---------------------------------------------------------------+---------------------+-----------------------------------------------------------------------------------------------+

4575
ncl_reference/psadilookup.dat
View File

File diff suppressed because it is too large Load Diff

215
ncl_reference/rcm2points.f
Unescape View File

@ -1,215 +0,0 @@
c -----------------------------------------------------------
C NCLFORTSTART
SUBROUTINE DRCM2POINTS(NGRD,NYI,NXI,YI,XI,FI,NXYO,YO,XO,FO
+ ,XMSG,OPT,NCRIT,KVAL,IER)
IMPLICIT NONE
INTEGER NGRD,NXI,NYI,NXYO,OPT,NCRIT,KVAL,IER
DOUBLE PRECISION XI(NXI,NYI),YI(NXI,NYI),FI(NXI,NYI,NGRD)
DOUBLE PRECISION XO(NXYO),YO(NXYO),FO(NXYO,NGRD),XMSG
C NCLEND
C This is written with GNU f77 acceptable extensions
c . This could be improved considerably with f90
c nomenclature:
c . nxi,nyi - lengths of xi,yi and dimensions of fi (must be >= 2)
c . xi - coordinates of fi (eg, lon [2D] )
c . yi - coordinates of fi (eg, lat [2D] )
c . fi - functional input values [2D]
c . nxyo - number of output points
c . xo - lon coordinates of fo (eg, lon [1D])
c . yo - lat coordinates of fo (eg, lat [1D])
c . fo - functional output values [interpolated]
c . xmsg - missing code
c . opt - 0/1 = inv distance, 2 = bilinear
c . ier - error code
c . =0; no error
c . =1; not enough points in input/output array
c . =2/3; xi or yi are not monotonically increasing
c . =4/5; xo or yo are not monotonically increasing
c
c local
INTEGER NG,NX,NY,NXY,NEXACT,IX,IY,M,N,NW,NER,K
DOUBLE PRECISION FW(2,2),W(2,2),SUMF,SUMW,CHKLAT(NYI),CHKLON(NXI)
DOUBLE PRECISION DGCDIST, WX, WY
DOUBLE PRECISION REARTH, DLAT, PI, RAD, DKM, DIST
c error checking
IER = 0
IF (NXI.LE.1 .OR. NYI.LE.1 .OR. NXYO.LE.0) THEN
IER = 1
RETURN
END IF
IF (IER.NE.0) RETURN
DO NY = 1,NYI
CHKLAT(NY) = YI(1,NY)
c c c print *,"chklat: ny=",ny," chklat=",chklat(ny)
END DO
CALL DMONOINC(CHKLAT,NYI,IER,NER)
IF (IER.NE.0) RETURN
DO NX = 1,NXI
CHKLON(NX) = XI(NX,1)
c c c print *,"chklon: nx=",nx," chklon=",chklon(nx)
END DO
CALL DMONOINC(CHKLAT,NYI,IER,NER)
IF (IER.NE.0) RETURN
C ORIGINAL (k = op, never implemented)
IF (KVAL.LE.0) THEN
K = 1
ELSE
K = KVAL
END IF
DO NG = 1,NGRD
DO NXY = 1,NXYO
FO(NXY,NG) = XMSG
END DO
END DO
c main loop [exact matches]
NEXACT = 0
DO NXY = 1,NXYO
DO IY = 1,NYI
DO IX = 1,NXI
IF (XO(NXY).EQ.XI(IX,IY) .AND.
+ YO(NXY).EQ.YI(IX,IY)) THEN
DO NG = 1,NGRD
FO(NXY,NG) = FI(IX,IY,NG)
NEXACT = NEXACT + 1
END DO
GO TO 10
END IF
END DO
END DO
10 CONTINUE
END DO
c c c print *, "nexact=",nexact
c main loop [interpolation]
DO NXY = 1,NXYO
DO IY = 1,NYI - K
DO IX = 1,NXI - K
IF (XO(NXY).GE.XI(IX,IY) .AND.
+ XO(NXY).LE.XI(IX+K,IY) .AND.
+ YO(NXY).GE.YI(IX,IY) .AND.
+ YO(NXY).LE.YI(IX,IY+K)) THEN
IF (ABS(OPT).EQ.2) THEN
WX = (XO(NXY)-XI(IX,IY))/
+ (XI(IX+K,IY)-XI(IX,IY))
WY = (YO(NXY)-YI(IX,IY))/
+ (YI(IX,IY+K)-YI(IX,IY))
W(1,1) = (1.D0-WX)*(1.D0-WY)
W(2,1) = WX*(1.D0-WY)
W(1,2) = (1.D0-WX)*WY
W(2,2) = WX*WY
ELSE
W(1,1) = (1.D0/DGCDIST(YO(NXY),XO(NXY),
+ YI(IX,IY),XI(IX,IY),2))**2
W(2,1) = (1.D0/DGCDIST(YO(NXY),XO(NXY),
+ YI(IX+K,IY),XI(IX+K,IY),2))**2
W(1,2) = (1.D0/DGCDIST(YO(NXY),XO(NXY),
+ YI(IX,IY+K),XI(IX,IY+K),2))**2
W(2,2) = (1.D0/DGCDIST(YO(NXY),XO(NXY),
+ YI(IX+K,IY+K),XI(IX+K,IY+K),2))**2
END IF
DO NG = 1,NGRD
IF (FO(NXY,NG).EQ.XMSG) THEN
FW(1,1) = FI(IX,IY,NG)
FW(2,1) = FI(IX+K,IY,NG)
FW(1,2) = FI(IX,IY+K,NG)
FW(2,2) = FI(IX+K,IY+K,NG)
NW = 0
SUMF = 0.0D0
SUMW = 0.0D0
DO N = 1,2
DO M = 1,2
IF (FW(M,N).NE.XMSG) THEN
SUMF = SUMF + FW(M,N)*W(M,N)
SUMW = SUMW + W(M,N)
NW = NW + 1
END IF
END DO
END DO
c nw >=3 arbitrary
IF (NW.GE.NCRIT .AND. SUMW.GT.0.D0) THEN
FO(NXY,NG) = SUMF/SUMW
END IF
END IF
END DO
GO TO 20
END IF
END DO
END DO
20 CONTINUE
END DO
C Are all the output points filled in? Check the 1st grid
C If so, return
DO NG = 1,NGRD
DO NXY = 1,NXYO
IF (FO(NXY,NG).EQ.XMSG) GO TO 30
END DO
END DO
RETURN
C only enter if some points are not interpolated to
C DLAT is arbitrary. It ould be made an option.
C DLAT is expressed in terms of degrees of latitude.
C DKM is DLAT in KILOMETERS
30 REARTH= 6371D0
DLAT = 5
PI = 4D0*ATAN(1.0D0)
RAD = PI/180D0
DKM = DLAT*(2D0*PI*REARTH)/360D0
C LOOP OVER EACH GRID ... INEFFICIENT
C THE RUB IS THAT SOME LEVELS COULD HAVE XMSG.
DO NG = 1,NGRD
DO NXY = 1,NXYO
IF(FO(NXY,NG).EQ.XMSG) THEN
C FIND ALL GRID POINTS WITHIN 'DKM' KILOMETERS OF PT
NW = 0
SUMF = 0.0D0
SUMW = 0.0D0
DO IY = 1,NYI
DO IX = 1,NXI
IF ((YI(IX,IY).GE.YO(NXY)-DLAT) .AND.
+ (YI(IX,IY).LE.YO(NXY)+DLAT)) THEN
DIST = DGCDIST(YO(NXY),XO(NXY)
+ ,YI(IX,IY),XI(IX,IY),2)
IF (DIST.LE.DKM .AND. DIST.GT.0.0D0 .AND.
+ FI(IX,IY,NG).NE.XMSG) THEN
DIST = 1.0D0/DIST**2
SUMF = SUMF + FI(IX,IY,NG)*DIST
SUMW = SUMW + DIST
NW = NW + 1
END IF
END IF
END DO
END DO
C C C IF (NW.GE.NCRIT .AND. SUMW.GT. 0.0D0) THEN
IF (SUMW.GT.0.0D0) THEN
FO(NXY,NG) = SUMF/SUMW
END IF
END IF
END DO
END DO
RETURN
END

376
ncl_reference/rcm2rgrid.f
Unescape View File

@ -1,376 +0,0 @@
C NCLFORTSTART
SUBROUTINE DRCM2RGRID(NGRD,NYI,NXI,YI,XI,FI,NYO,YO,NXO,XO,FO
+ ,XMSG,NCRIT,OPT,IER)
IMPLICIT NONE
INTEGER NGRD,NXI,NYI,NXO,NYO,NCRIT,OPT,IER
DOUBLE PRECISION XI(NXI,NYI),YI(NXI,NYI),FI(NXI,NYI,NGRD)
DOUBLE PRECISION XO(NXO),YO(NYO),FO(NXO,NYO,NGRD),XMSG
C NCLEND
C This is written with GNU f77 acceptable extensions
c . This could be improved considerably with f90
c NCL: fo = rcm2rgrid (lat2d,lon2d,fi, lat, lon iopt)
c yi xi fi yo xo
c
c fo is the same size xo, yo and same type as "fi"
c xmsg = fi@_FillValue
c opt unused option
c
c The NCL wrapper should allow for multiple datasets
c so the user need only make one call to the function.
c perform 2D interpolation allowing for missing data: nothing fancy
c nomenclature:
c . nxi,nyi - lengths of xi,yi and dimensions of fi (must be >= 2)
c . xi - coordinates of fi (eg, lon [2D] )
c . yi - coordinates of fi (eg, lat [2D] )
c . fi - functional input values [2D]
c . nxo,nyo - lengths of xo,yo and dimensions of fo (must be >= 1)
c . xo - coordinates of fo (eg, lon [1D])
c . must be monotonically increasing
c . yo - coordinates of fo (eg, lat [1D])
c . must be monotonically increasing
c . fo - functional output values [interpolated]
c . xmsg - missing code
c . opt - unused
c . ier - error code
c . =0; no error
c . =1; not enough points in input/output array
c . =2/3; xi or yi are not monotonically increasing
c . =4/5; xo or yo are not monotonically increasing
c
c local
INTEGER NG, NX,NY,NEXACT,IX,IY,M,N,NW,NER,K,NCRT
INTEGER MFLAG, MPTCRT, MKNT
DOUBLE PRECISION FW(2,2),W(2,2),SUMF,SUMW,CHKLAT(NYI),CHKLON(NXI)
DOUBLE PRECISION EPS
DOUBLE PRECISION DGCDIST
c error checking
IER = 0
IF (NXI.LE.1 .OR. NYI.LE.1 .OR. NXO.LE.1 .OR. NYO.LE.1) THEN
IER = 1
RETURN
END IF
IF (IER.NE.0) RETURN
CALL DMONOINC(YO,NYO,IER,NER)
IF (IER.NE.0) RETURN
CALL DMONOINC(XO,NXO,IER,NER)
IF (IER.NE.0) RETURN
DO NY = 1,NYI
CHKLAT(NY) = YI(1,NY)
c c c print *,"chklat: ny=",ny," chklat=",chklat(ny)
END DO
CALL DMONOINC(CHKLAT,NYI,IER,NER)
IF (IER.NE.0) RETURN
DO NX = 1,NXI
CHKLON(NX) = XI(NX,1)
c c c print *,"chklon: nx=",nx," chklon=",chklon(nx)
END DO
CALL DMONOINC(CHKLAT,NYI,IER,NER)
IF (IER.NE.0) RETURN
K = 2
c c c k = opt
IF (NCRIT.LE.1) THEN
NCRT = 1
ELSE
NCRT = MIN(4,NCRIT)
END IF
c initialize to xmsg
DO NG=1,NGRD
DO NY = 1,NYO
DO NX = 1,NXO
FO(NX,NY,NG) = XMSG
END DO
END DO
END DO
c main loop [exact matches]
c people want bit-for-bit match
EPS = 1.D-04
NEXACT = 0
DO NY = 1,NYO
DO NX = 1,NXO
DO IY = 1,NYI
DO IX = 1,NXI
IF (XO(NX).GE.(XI(IX,IY)-EPS) .AND.
+ XO(NX).LE.(XI(IX,IY)+EPS) .AND.
+ YO(NY).GE.(YI(IX,IY)-EPS) .AND.
+ YO(NY).LE.(YI(IX,IY)+EPS) ) THEN
DO NG=1,NGRD
FO(NX,NY,NG) = FI(IX,IY,NG)
NEXACT = NEXACT + 1
END DO
GO TO 10
END IF
END DO
END DO
10 CONTINUE
END DO
END DO
c c c print *, "nexact=",nexact
c main loop [interpolation]
DO NY = 1,NYO
DO NX = 1,NXO
DO IY = 1,NYI-K
DO IX = 1,NXI-K
IF (XO(NX).GE.XI(IX,IY) .AND.
+ XO(NX).LE.XI(IX+K,IY) .AND.
+ YO(NY).GE.YI(IX,IY) .AND.
+ YO(NY).LE.YI(IX,IY+K)) THEN
W(1,1) = (1.D0/DGCDIST(YO(NY),XO(NX),
+ YI(IX,IY),XI(IX,IY),2))**2
W(2,1) = (1.D0/DGCDIST(YO(NY),XO(NX),
+ YI(IX+K,IY),XI(IX+K,IY),2))**2
W(1,2) = (1.D0/DGCDIST(YO(NY),XO(NX),
+ YI(IX,IY+K),XI(IX,IY+K),2))**2
W(2,2) = (1.D0/DGCDIST(YO(NY),XO(NX),
+ YI(IX+K,IY+K),XI(IX+K,IY+K),2))**2
DO NG=1,NGRD
IF (FO(NX,NY,NG).EQ.XMSG) THEN
FW(1,1) = FI(IX,IY,NG)
FW(2,1) = FI(IX+K,IY,NG)
FW(1,2) = FI(IX,IY+K,NG)
FW(2,2) = FI(IX+K,IY+K,NG)
NW = 0
SUMF = 0.0D0
SUMW = 0.0D0
DO N = 1,2
DO M = 1,2
IF (FW(M,N).NE.XMSG) THEN
SUMF = SUMF + FW(M,N)*W(M,N)
SUMW = SUMW + W(M,N)
NW = NW + 1
END IF
END DO
END DO
c nw >=3 arbitrary
c c c IF (NW.GE.3 .AND. SUMW.GT.0.D0) THEN
c nw =1 nearest neighbor
IF (NW.GE.NCRT .AND. SUMW.GT.0.D0) THEN
FO(NX,NY,NG) = SUMF/SUMW
END IF
END IF
END DO
GO TO 20
END IF
END DO
END DO
20 CONTINUE
END DO
END DO
C Since the RCM grid is curvilinear the above algorithm may not work
C . for all of the locations on regular grid. Fill via linear interp.
MKNT = 0
MFLAG = 0
MPTCRT = 2
DO NG=1,NGRD
DO NY=1,NYO
DO NX=1,NXO
IF (FO(NX,NY,NG).EQ.XMSG) THEN
CALL DLINMSG(FO(1,NY,NG),NXO,XMSG,MFLAG,MPTCRT)
MKNT = MKNT + 1
END IF
END DO
END DO
END DO
C C C PRINT *,"MKNT=",MKNT
RETURN
END
c -----------------------------------------------------------
C NCLFORTSTART
SUBROUTINE DRGRID2RCM(NGRD,NYI,NXI,YI,XI,FI,NYO,NXO,YO,XO,FO
+ ,XMSG,NCRIT,OPT,IER)
IMPLICIT NONE
INTEGER NGRD,NXI,NYI,NXO,NYO,OPT,NCRIT,IER
DOUBLE PRECISION XI(NXI),YI(NYI),FI(NXI,NYI,NGRD)
DOUBLE PRECISION XO(NXO,NYO),YO(NXO,NYO),FO(NXO,NYO,NGRD),XMSG
C NCLEND
C This is written with GNU f77 acceptable extensions
c . This could be improved considerably with f90
c fo is the same size xo, yo and same type as "fi"
c xmsg = fi@_FillValue
c opt unused option
c
c The NCL wrapper should allow for multiple datasets
c so the user need only make one call to the function.
c perform 2D interpolation allowing for missing data: nothing fancy
c nomenclature:
c . nxi,nyi - lengths of xi,yi and dimensions of fi (must be >= 2)
c . xi - coordinates of fi (eg, lon [1D])
c . yi - coordinates of fi (eg, lat [1D])
c . fi - functional input values [2D]
c . nxo,nyo - lengths of xo,yo and dimensions of fo (must be >= 1)
c . xo - coordinates of fo (eg, lon [2D])
c . must be monotonically increasing
c . yo - coordinates of fo (eg, lat [2D])
c . must be monotonically increasing
c . fo - functional output values [interpolated]
c . xmsg - missing code
c . opt - unused
c . ier - error code
c . =0; no error
c . =1; not enough points in input/output array
c . =2/3; xi or yi are not monotonically increasing
c . =4/5; xo or yo are not monotonically increasing
c
c local
INTEGER NG,NX,NY,NEXACT,IX,IY,M,N,NW,NER,K
DOUBLE PRECISION FW(2,2),W(2,2),SUMF,SUMW,EPS
DOUBLE PRECISION DGCDIST
c in-line functions (bilinear interp)
DOUBLE PRECISION Z1,Z2,Z3,Z4,SLOPE,SLPX,SLPY,FLI,FBLI
FLI(Z1,Z2,SLOPE) = Z1 + SLOPE* (Z2-Z1)
FBLI(Z1,Z2,Z3,Z4,SLPX,SLPY) = FLI(Z1,Z2,SLPX) +
+ SLPY* (FLI(Z3,Z4,SLPX)-
+ FLI(Z1,Z2,SLPX))
c error checking
IER = 0
IF (NXI.LE.1 .OR. NYI.LE.1 .OR. NXO.LE.1 .OR. NYO.LE.1) THEN
IER = 1
RETURN
END IF
IF (IER.NE.0) RETURN
CALL DMONOINC(YI,NYI,IER,NER)
IF (IER.NE.0) RETURN
CALL DMONOINC(XI,NXI,IER,NER)
IF (IER.NE.0) RETURN
c Init to missing
DO NG = 1,NGRD
DO NY = 1,NYO
DO NX = 1,NXO
FO(NX,NY,NG) = XMSG
END DO
END DO
END DO
c main loop [exact matches]
EPS = 1.D-03
NEXACT = 0
DO NY = 1,NYO
DO NX = 1,NXO
DO IY = 1,NYI
DO IX = 1,NXI
IF (XO(NX,NY).GE.(XI(IX)-EPS) .AND.
+ XO(NX,NY).LE.(XI(IX)+EPS) .AND.
+ YO(NX,NY).GE.(YI(IY)-EPS) .AND.
+ YO(NX,NY).LE.(YI(IY)+EPS) ) THEN
DO NG=1,NGRD
FO(NX,NY,NG) = FI(IX,IY,NG)
NEXACT = NEXACT + 1
END DO
GO TO 10
END IF
END DO
END DO
10 CONTINUE
END DO
END DO
c c c print *, "nexact=",nexact
K = 1
c c c k = opt
c main loop [interpolation]
DO NY = 1,NYO
DO NX = 1,NXO
DO IY = 1,NYI - K
DO IX = 1,NXI - K
IF (XO(NX,NY).GE.XI(IX) .AND.
+ XO(NX,NY).LT.XI(IX+K) .AND.
+ YO(NX,NY).GE.YI(IY) .AND.
+ YO(NX,NY).LT.YI(IY+K)) THEN
DO NG = 1,NGRD
IF (FO(NX,NY,NG).EQ.XMSG) THEN
IF (FI(IX,IY,NG).NE.XMSG .AND.
+ FI(IX+K,IY,NG).NE.XMSG .AND.
+ FI(IX,IY+K,NG).NE.XMSG .AND.
+ FI(IX+K,IY+K,NG).NE.XMSG) THEN
FO(NX,NY,NG) =FBLI(FI(IX,IY,NG),FI(IX+K,IY,NG),
+ FI(IX,IY+K,NG),FI(IX+K,IY+K,NG),
+ (XO(NX,NY)-XI(IX))/
+ (XI(IX+K)-XI(IX)),
+ (YO(NX,NY)-YI(IY))/
+ (YI(IY+K)-YI(IY)))
ELSE
c OVERKILL
FW(1,1) = FI(IX,IY,NG)
FW(2,1) = FI(IX+K,IY,NG)
FW(1,2) = FI(IX,IY+K,NG)
FW(2,2) = FI(IX+K,IY+K,NG)
W(1,1) = (1.D0/DGCDIST(YO(NX,NY),XO(NX,NY)
+ ,YI(IY),XI(IX),2))**2
W(2,1) = (1.D0/DGCDIST(YO(NX,NY),XO(NX,NY)
+ ,YI(IY),XI(IX+K),2))**2
W(1,2) = (1.D0/DGCDIST(YO(NX,NY),XO(NX,NY)
+ ,YI(IY+K),XI(IX),2))**2
W(2,2) = (1.D0/DGCDIST(YO(NX,NY),XO(NX,NY)
+ ,YI(IY+K),XI(IX+K),2))**2
NW = 0
SUMF = 0.0D0
SUMW = 0.0D0
DO N = 1,2
DO M = 1,2
IF (FW(M,N).NE.XMSG) THEN
SUMF = SUMF + FW(M,N)*W(M,N)
SUMW = SUMW + W(M,N)
NW = NW + 1
END IF
END DO
END DO
c nw >=3 arbitrary
c c c IF (NCRIT.GE.3 .AND. SUMW.GT.0.D0) THEN
c nw =1 nearest neighbor
IF (NCRIT.GE.1 .AND. SUMW.GT.0.D0) THEN
FO(NX,NY,NG) = SUMF/SUMW
END IF
END IF
END IF
END DO
GO TO 20
END IF
END DO
END DO
20 CONTINUE
END DO
END DO
RETURN
END

832
ncl_reference/rcmW.c
Unescape View File

@ -1,832 +0,0 @@
#include <stdio.h>
#include "wrapper.h"
extern void NGCALLF(drcm2rgrid,DRCM2RGRID)(int *,int *,int *,double *,double *,
double *,int *,double *,int*,
double *,double *,double*,
int *,int *,int *);
extern void NGCALLF(drgrid2rcm,DRGRID2RCM)(int *,int *,int *,double *,double *,
double *,int *,int *,double *,
double *,double *,double*,
int *,int *,int *);
extern void NGCALLF(drcm2points,DRCM2POINTS)(int *,int *,int *,double *,
double *,double *,int *,double *,
double *,double *,double*,
int *,int *,int *,int*);
NhlErrorTypes rcm2rgrid_W( void )
{
/*
* Input variables
*/
void *lat2d, *lon2d, *fi, *lat1d, *lon1d, *opt;
double *tmp_lat2d, *tmp_lon2d, *tmp_lat1d, *tmp_lon1d, *tmp_fi;
int tmp_opt, tmp_ncrit;
ng_size_t dsizes_lat2d[2];
ng_size_t dsizes_lon2d[2];
ng_size_t dsizes_lat1d[2];
ng_size_t dsizes_lon1d[1];
int ndims_fi;
ng_size_t size_fi;
ng_size_t dsizes_fi[NCL_MAX_DIMENSIONS];
int has_missing_fi;
NclScalar missing_fi, missing_dfi, missing_rfi;
NclBasicDataTypes type_lat2d, type_lon2d, type_lat1d, type_lon1d;
NclBasicDataTypes type_fi, type_opt;
/*
* Output variables.
*/
void *fo;
double *tmp_fo;
ng_size_t *dsizes_fo;
NclBasicDataTypes type_fo;
NclScalar missing_fo;
/*
* Other variables
*/
ng_size_t nlon2d, nlat2d, nfi, nlat1d, nlon1d, nfo, ngrid, size_fo;
ng_size_t i;
int ier, ret;
int inlon2d, inlat2d, ingrid, inlon1d, inlat1d;
/*
* Retrieve parameters
*
* Note that any of the pointer parameters can be set to NULL,
* which implies you don't care about its value.
*/
lat2d = (void*)NclGetArgValue(
0,
6,
NULL,
dsizes_lat2d,
NULL,
NULL,
&type_lat2d,
DONT_CARE);
lon2d = (void*)NclGetArgValue(
1,
6,
NULL,
dsizes_lon2d,
NULL,
NULL,
&type_lon2d,
DONT_CARE);
fi = (void*)NclGetArgValue(
2,
6,
&ndims_fi,
dsizes_fi,
&missing_fi,
&has_missing_fi,
&type_fi,
DONT_CARE);
lat1d = (void*)NclGetArgValue(
3,
6,
NULL,
dsizes_lat1d,
NULL,
NULL,
&type_lat1d,
DONT_CARE);
lon1d = (void*)NclGetArgValue(
4,
6,
NULL,
dsizes_lon1d,
NULL,
NULL,
&type_lon1d,
DONT_CARE);
opt = (void*)NclGetArgValue(
5,
6,
NULL,
NULL,
NULL,
NULL,
&type_opt,
DONT_CARE);
/*
* Check the input lat/lon arrays. They must be the same size, and larger
* than one element.
*/
if(dsizes_lat2d[0] != dsizes_lon2d[0] ||
dsizes_lat2d[1] != dsizes_lon2d[1]) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rcm2rgrid: The input lat/lon grids must be the same size");
return(NhlFATAL);
}
nlat2d = dsizes_lat2d[0];
nlon2d = dsizes_lat2d[1]; /* same as dsizes_lon2d[1] */
nlat1d = dsizes_lat1d[0];
nlon1d = dsizes_lon1d[0];
if(nlon2d <= 1 || nlat2d <= 1 || nlat1d <= 1 || nlon1d <= 1) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rcm2rgrid: The input/output lat/lon grids must have at least 2 elements");
return(NhlFATAL);
}
/*
* Compute the total number of elements in our arrays.
*/
nfi = nlon2d * nlat2d;
nfo = nlat1d * nlon1d;
/*
* Check dimensions of fi.
*/
if(ndims_fi < 2) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rcm2rgrid: fi must be at least two dimensions");
return(NhlFATAL);
}
if(dsizes_fi[ndims_fi-2] != nlat2d || dsizes_fi[ndims_fi-1] != nlon2d) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rcm2rgrid: The rightmost dimensions of fi must be nlat2d x nlon2d, where nlat2d and nlon2d are the dimensions of the lat2d/lon2d arrays");
return(NhlFATAL);
}
/*
* Compute the total size of the input/output arrays.
*/
ngrid = 1;
for( i = 0; i < ndims_fi-2; i++ ) ngrid *= dsizes_fi[i];
size_fi = ngrid * nfi;
size_fo = ngrid * nfo;
/*
* Test input dimension sizes.
*/
if((nlon2d > INT_MAX) || (nlat2d > INT_MAX) || (ngrid > INT_MAX) ||
(nlon1d > INT_MAX) || (nlat1d > INT_MAX)) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rcm2rgrid: one or more input dimension sizes is greater than INT_MAX");
return(NhlFATAL);
}
inlon2d = (int) nlon2d;
inlat2d = (int) nlat2d;
ingrid = (int) ngrid;
inlon1d = (int) nlon1d;
inlat1d = (int) nlat1d;
/*
* Coerce missing values.
*/
coerce_missing(type_fi,has_missing_fi,&missing_fi,&missing_dfi,
&missing_rfi);
/*
* Allocate space for output array.
*/
if(type_fi == NCL_double) {
fo = (void*)calloc(size_fo,sizeof(double));
tmp_fo = &((double*)fo)[0];
type_fo = NCL_double;
missing_fo.doubleval = missing_dfi.doubleval;
}
else {
fo = (void*)calloc(size_fo,sizeof(float));
tmp_fo = (double*)calloc(size_fo,sizeof(double));
type_fo = NCL_float;
missing_fo.floatval = missing_rfi.floatval;
if(tmp_fo == NULL) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rcm2rgrid: Unable to allocate memory for temporary array");
return(NhlFATAL);
}
}
dsizes_fo = (ng_size_t*)calloc(ndims_fi,sizeof(ng_size_t));
if(fo == NULL || dsizes_fo == NULL) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rcm2rgrid: Unable to allocate memory for output array");
return(NhlFATAL);
}
for(i = 0; i < ndims_fi-2; i++) dsizes_fo[i] = dsizes_fi[i];
dsizes_fo[ndims_fi-2] = nlat1d;
dsizes_fo[ndims_fi-1] = nlon1d;
/*
* Coerce input arrays to double if necessary.
*/
tmp_lat2d = coerce_input_double(lat2d,type_lat2d,nfi,0,NULL,NULL);
tmp_lon2d = coerce_input_double(lon2d,type_lon2d,nfi,0,NULL,NULL);
tmp_lat1d = coerce_input_double(lat1d,type_lat1d,nlat1d,0,NULL,NULL);
tmp_lon1d = coerce_input_double(lon1d,type_lon1d,nlon1d,0,NULL,NULL);
tmp_fi = coerce_input_double(fi,type_fi,size_fi,has_missing_fi,
&missing_fi,&missing_dfi);
if(tmp_lat2d == NULL || tmp_lon2d == NULL ||
tmp_lat1d == NULL || tmp_lon1d == NULL || tmp_fi == NULL) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rcm2rgrid: Unable to coerce input lat/lon arrays to double precision");
return(NhlFATAL);
}
/*
* Force opt to zero and ncrit to 1, since they are not used yet.
*/
tmp_opt = 0;
tmp_ncrit = 1;
NGCALLF(drcm2rgrid,DRCM2RGRID)(&ingrid,&inlat2d,&inlon2d,tmp_lat2d,tmp_lon2d,
tmp_fi,&inlat1d,tmp_lat1d,&inlon1d,
tmp_lon1d,tmp_fo,&missing_dfi.doubleval,
&tmp_ncrit,&tmp_opt,&ier);
if(ier) {
if(ier == 1) {
NhlPError(NhlWARNING,NhlEUNKNOWN,"rcm2rgrid: not enough points in input/output array");
}
if(2 <= ier && ier <= 5) {
NhlPError(NhlWARNING,NhlEUNKNOWN,"rcm2rgrid: lat2d, lon2d, lat1d, lon1d must be monotonically increasing");
}
set_subset_output_missing(fo,0,type_fo,size_fo,missing_dfi.doubleval);
}
else {
if(type_fo != NCL_double) {
coerce_output_float_only(fo,tmp_fo,size_fo,0);
}
}
/*
* Free temp arrays.
*/
if(type_lat2d != NCL_double) NclFree(tmp_lat2d);
if(type_lon2d != NCL_double) NclFree(tmp_lon2d);
if(type_lat1d != NCL_double) NclFree(tmp_lat1d);
if(type_lon1d != NCL_double) NclFree(tmp_lon1d);
if(type_fi != NCL_double) NclFree(tmp_fi);
if(type_fo != NCL_double) NclFree(tmp_fo);
/*
* Return.
*/
ret = NclReturnValue(fo,ndims_fi,dsizes_fo,&missing_fo,type_fo,0);
NclFree(dsizes_fo);
return(ret);
}
NhlErrorTypes rgrid2rcm_W( void )
{
/*
* Input variables
*/
void *lat2d, *lon2d, *fi, *lat1d, *lon1d, *opt;
double *tmp_lat2d, *tmp_lon2d, *tmp_lat1d, *tmp_lon1d, *tmp_fi;
int tmp_opt, tmp_ncrit;
ng_size_t dsizes_lat2d[2];
ng_size_t dsizes_lon2d[2];
ng_size_t dsizes_lat1d[2];
ng_size_t dsizes_lon1d[1];
int ndims_fi;
ng_size_t dsizes_fi[NCL_MAX_DIMENSIONS];
int has_missing_fi;
NclScalar missing_fi, missing_dfi, missing_rfi;
NclBasicDataTypes type_lat2d, type_lon2d, type_lat1d, type_lon1d;
NclBasicDataTypes type_fi, type_opt;
/*
* Output variables.
*/
void *fo;
double *tmp_fo;
ng_size_t *dsizes_fo;
NclBasicDataTypes type_fo;
NclScalar missing_fo;
/*
* Other variables
*/
ng_size_t nlon2d, nlat2d, nfi, nlat1d, nlon1d, nfo, ngrid, size_fi, size_fo;
ng_size_t i;
int ier, ret;
int inlon2d, inlat2d, ingrid, inlon1d, inlat1d;
/*
* Retrieve parameters
*
* Note that any of the pointer parameters can be set to NULL,
* which implies you don't care about its value.
*/
lat1d = (void*)NclGetArgValue(
0,
6,
NULL,
dsizes_lat1d,
NULL,
NULL,
&type_lat1d,
DONT_CARE);
lon1d = (void*)NclGetArgValue(
1,
6,
NULL,
dsizes_lon1d,
NULL,
NULL,
&type_lon1d,
DONT_CARE);
fi = (void*)NclGetArgValue(
2,
6,
&ndims_fi,
dsizes_fi,
&missing_fi,
&has_missing_fi,
&type_fi,
DONT_CARE);
lat2d = (void*)NclGetArgValue(
3,
6,
NULL,
dsizes_lat2d,
NULL,
NULL,
&type_lat2d,
DONT_CARE);
lon2d = (void*)NclGetArgValue(
4,
6,
NULL,
dsizes_lon2d,
NULL,
NULL,
&type_lon2d,
DONT_CARE);
opt = (void*)NclGetArgValue(
5,
6,
NULL,
NULL,
NULL,
NULL,
&type_opt,
DONT_CARE);
/*
* Check the output lat/lon arrays. They must be the same size, and larger
* than one element.
*/
if(dsizes_lat2d[0] != dsizes_lon2d[0] ||
dsizes_lat2d[1] != dsizes_lon2d[1]) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rgrid2rcm: The output lat/lon grids must be the same size");
return(NhlFATAL);
}
nlat2d = dsizes_lat2d[0];
nlon2d = dsizes_lat2d[1]; /* same as dsizes_lon2d[1] */
nlat1d = dsizes_lat1d[0];
nlon1d = dsizes_lon1d[0];
if(nlon2d <= 1 || nlat2d <= 1 || nlat1d <= 1 || nlon1d <= 1) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rgrid2rcm: The input/output lat/lon grids must have at least 2 elements");
return(NhlFATAL);
}
/*
* Compute the total number of elements in our arrays.
*/
nfi = nlat1d * nlon1d;
nfo = nlon2d * nlat2d;
/*
* Check dimensions of fi.
*/
if(ndims_fi < 2) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rgrid2rcm: fi must be at least two dimensions");
return(NhlFATAL);
}
if(dsizes_fi[ndims_fi-2] != nlat1d || dsizes_fi[ndims_fi-1] != nlon1d) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rgrid2rcm: The rightmost dimensions of fi must be nlat1d x nlon1d, where nlat1d and nlon1d are the dimensions of the lat1d/lon1d arrays");
return(NhlFATAL);
}
/*
* Compute the total size of the input/output arrays.
*/
ngrid = 1;
for( i = 0; i < ndims_fi-2; i++ ) ngrid *= dsizes_fi[i];
size_fi = ngrid * nfi;
size_fo = ngrid * nfo;
/*
* Test input dimension sizes.
*/
if((nlon2d > INT_MAX) || (nlat2d > INT_MAX) || (ngrid > INT_MAX) ||
(nlon1d > INT_MAX) || (nlat1d > INT_MAX)) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rgrid2rcm: one or more input dimension sizes is greater than INT_MAX");
return(NhlFATAL);
}
inlon2d = (int) nlon2d;
inlat2d = (int) nlat2d;
ingrid = (int) ngrid;
inlon1d = (int) nlon1d;
inlat1d = (int) nlat1d;
/*
* Coerce missing values.
*/
coerce_missing(type_fi,has_missing_fi,&missing_fi,&missing_dfi,
&missing_rfi);
/*
* Allocate space for output array.
*/
if(type_fi == NCL_double) {
fo = (void*)calloc(size_fo,sizeof(double));
tmp_fo = &((double*)fo)[0];
type_fo = NCL_double;
missing_fo.doubleval = missing_dfi.doubleval;
}
else {
tmp_fo = (double*)calloc(size_fo,sizeof(double));
fo = (void*)calloc(size_fo,sizeof(float));
type_fo = NCL_float;
missing_fo.floatval = missing_rfi.floatval;
if(tmp_fo == NULL) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rgrid2rcm: Unable to allocate memory for temporary arrays");
return(NhlFATAL);
}
}
dsizes_fo = (ng_size_t*)calloc(ndims_fi,sizeof(ng_size_t));
if(fo == NULL || dsizes_fo == NULL) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rgrid2rcm: Unable to allocate memory for output array");
return(NhlFATAL);
}
for(i = 0; i < ndims_fi-2; i++) dsizes_fo[i] = dsizes_fi[i];
dsizes_fo[ndims_fi-2] = nlat2d;
dsizes_fo[ndims_fi-1] = nlon2d;
/*
* Coerce input arrays to double if necessary.
*/
tmp_lat2d = coerce_input_double(lat2d,type_lat2d,nfo,0,NULL,NULL);
tmp_lon2d = coerce_input_double(lon2d,type_lon2d,nfo,0,NULL,NULL);
tmp_lat1d = coerce_input_double(lat1d,type_lat1d,nlat1d,0,NULL,NULL);
tmp_lon1d = coerce_input_double(lon1d,type_lon1d,nlon1d,0,NULL,NULL);
tmp_fi = coerce_input_double(fi,type_fi,size_fi,has_missing_fi,
&missing_fi,&missing_dfi);
if(tmp_lat2d == NULL || tmp_lon2d == NULL ||
tmp_lat1d == NULL || tmp_lon1d == NULL || tmp_fi == NULL) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rgrid2rcm: Unable to coerce input lat/lon arrays to double precision");
return(NhlFATAL);
}
/*
* Force opt to zero and ncrit to 1, since they are not used yet.
*/
tmp_opt = 0;
tmp_ncrit = 1;
NGCALLF(drgrid2rcm,DRGRID2RCM)(&ingrid,&inlat1d,&inlon1d,tmp_lat1d,tmp_lon1d,
tmp_fi,&inlat2d,&inlon2d,tmp_lat2d,
tmp_lon2d,tmp_fo,&missing_dfi.doubleval,
&tmp_ncrit,&tmp_opt,&ier);
if(ier) {
if(ier == 1) {
NhlPError(NhlWARNING,NhlEUNKNOWN,"rgrid2rcm: not enough points in input/output array");
}
if(2 <= ier && ier <= 5) {
NhlPError(NhlWARNING,NhlEUNKNOWN,"rgrid2rcm: lat2d, lon2d, lat1d, lon1d must be monotonically increasing");
}
set_subset_output_missing(fo,0,type_fo,size_fo,missing_dfi.doubleval);
}
else {
if(type_fo != NCL_double) {
coerce_output_float_only(fo,tmp_fo,size_fo,0);
}
}
/*
* Free temp arrays.
*/
if(type_lat2d != NCL_double) NclFree(tmp_lat2d);
if(type_lon2d != NCL_double) NclFree(tmp_lon2d);
if(type_lat1d != NCL_double) NclFree(tmp_lat1d);
if(type_lon1d != NCL_double) NclFree(tmp_lon1d);
if(type_fi != NCL_double) NclFree(tmp_fi);
if(type_fo != NCL_double) NclFree(tmp_fo);
ret = NclReturnValue(fo,ndims_fi,dsizes_fo,&missing_fo,type_fo,0);
NclFree(dsizes_fo);
return(ret);
}
NhlErrorTypes rcm2points_W( void )
{
/*
* Input variables
*/
void *lat2d, *lon2d, *fi, *lat1d, *lon1d;
double *tmp_lat2d, *tmp_lon2d, *tmp_lat1d, *tmp_lon1d, *tmp_fi;
int *opt, tmp_ncrit;
ng_size_t dsizes_lat2d[2];
ng_size_t dsizes_lon2d[2];
ng_size_t dsizes_lat1d[2];
ng_size_t dsizes_lon1d[1];
int ndims_fi;
ng_size_t dsizes_fi[NCL_MAX_DIMENSIONS];
int has_missing_fi;
NclScalar missing_fi, missing_dfi, missing_rfi;
NclBasicDataTypes type_lat2d, type_lon2d, type_lat1d, type_lon1d;
NclBasicDataTypes type_fi;
/*
* Variables for retrieving attributes from "opt".
*/
NclAttList *attr_list;
NclAtt attr_obj;
NclStackEntry stack_entry;
logical set_search_width;
/*
* Output variables.
*/
void *fo;
double *tmp_fo;
int ndims_fo;
ng_size_t *dsizes_fo;
NclBasicDataTypes type_fo;
NclScalar missing_fo;
/*
* Other variables
*/
ng_size_t nlon2d, nlat2d, nfi, nlat1d, nfo, ngrid, size_fi, size_fo;
ng_size_t i;
int search_width, ier, ret;
int inlon2d, inlat2d, ingrid, inlat1d;
/*
* Retrieve parameters
*
* Note that any of the pointer parameters can be set to NULL,
* which implies you don't care about its value.
*/
lat2d = (void*)NclGetArgValue(
0,
6,
NULL,
dsizes_lat2d,
NULL,
NULL,
&type_lat2d,
DONT_CARE);
lon2d = (void*)NclGetArgValue(
1,
6,
NULL,
dsizes_lon2d,
NULL,
NULL,
&type_lon2d,
DONT_CARE);
fi = (void*)NclGetArgValue(
2,
6,
&ndims_fi,
dsizes_fi,
&missing_fi,
&has_missing_fi,
&type_fi,
DONT_CARE);
lat1d = (void*)NclGetArgValue(
3,
6,
NULL,
dsizes_lat1d,
NULL,
NULL,
&type_lat1d,
DONT_CARE);
lon1d = (void*)NclGetArgValue(
4,
6,
NULL,
dsizes_lon1d,
NULL,
NULL,
&type_lon1d,
DONT_CARE);
opt = (int*)NclGetArgValue(
5,
6,
NULL,
NULL,
NULL,
NULL,
NULL,
DONT_CARE);
/*
* Check the input lat/lon arrays. They must be the same size, and larger
* than one element.
*/
if(dsizes_lat2d[0] != dsizes_lon2d[0] ||
dsizes_lat2d[1] != dsizes_lon2d[1]) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rcm2points: The input lat/lon grids must be the same size");
return(NhlFATAL);
}
nlat2d = dsizes_lat2d[0];
nlon2d = dsizes_lat2d[1]; /* same as dsizes_lon2d[1] */
nlat1d = dsizes_lat1d[0];
if(dsizes_lon1d[0] != nlat1d) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rcm2points: The output lat/lon arrays must be the same length");
return(NhlFATAL);
}
if(nlon2d < 2 || nlat2d < 2 || nlat1d < 1) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rcm2points: The input lat/lon grids must have at least 2 elements, and the output lat/lon arrays 1 element");
return(NhlFATAL);
}
/*
* Compute the total number of elements in our arrays.
*/
nfi = nlon2d * nlat2d;
nfo = nlat1d;
/*
* Check dimensions of fi.
*/
if(ndims_fi < 2) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rcm2points: fi must be at least two dimensions");
return(NhlFATAL);
}
if(dsizes_fi[ndims_fi-2] != nlat2d || dsizes_fi[ndims_fi-1] != nlon2d) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rcm2points: The rightmost dimensions of fi must be nlat2d x nlon2d, where nlat2d and nlon2d are the dimensions of the lat2d/lon2d arrays");
return(NhlFATAL);
}
/*
* Compute the sizes of the input/output arrays.
*/
ngrid = 1;
for( i = 0; i < ndims_fi-2; i++ ) ngrid *= dsizes_fi[i];
size_fi = ngrid * nfi;
size_fo = ngrid * nfo;
/*
* Test input dimension sizes.
*/
if((nlon2d > INT_MAX) || (nlat2d > INT_MAX) || (ngrid > INT_MAX) || (nlat1d > INT_MAX)) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rcm2points: one or more input dimension sizes is greater than INT_MAX");
return(NhlFATAL);
}
inlon2d = (int) nlon2d;
inlat2d = (int) nlat2d;
ingrid = (int) ngrid;
inlat1d = (int) nlat1d;
/*
* Coerce missing values.
*/
coerce_missing(type_fi,has_missing_fi,&missing_fi,&missing_dfi,
&missing_rfi);
/*
* Allocate space for output array.
*/
if(type_fi == NCL_double) {
fo = (void*)calloc(size_fo,sizeof(double));
tmp_fo = &((double*)fo)[0];
type_fo = NCL_double;
missing_fo.doubleval = missing_dfi.doubleval;
}
else {
fo = (void*)calloc(size_fo,sizeof(float));
tmp_fo = (double*)calloc(size_fo,sizeof(double));
type_fo = NCL_float;
missing_fo.floatval = missing_rfi.floatval;
if(tmp_fo == NULL) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rcm2points: Unable to allocate memory for temporary array");
return(NhlFATAL);
}
}
ndims_fo = ndims_fi-1;
dsizes_fo = (ng_size_t*)calloc(ndims_fo,sizeof(ng_size_t));
if(fo == NULL || dsizes_fo == NULL) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rcm2points: Unable to allocate memory for output array");
return(NhlFATAL);
}
for(i = 0; i < ndims_fi-2; i++) dsizes_fo[i] = dsizes_fi[i];
dsizes_fo[ndims_fi-2] = nlat1d;
/*
* Coerce input arrays to double if necessary.
*/
tmp_lat2d = coerce_input_double(lat2d,type_lat2d,nfi,0,NULL,NULL);
tmp_lon2d = coerce_input_double(lon2d,type_lon2d,nfi,0,NULL,NULL);
tmp_lat1d = coerce_input_double(lat1d,type_lat1d,nlat1d,0,NULL,NULL);
tmp_lon1d = coerce_input_double(lon1d,type_lon1d,nlat1d,0,NULL,NULL);
tmp_fi = coerce_input_double(fi,type_fi,size_fi,has_missing_fi,
&missing_fi,&missing_dfi);
if(tmp_lat2d == NULL || tmp_lon2d == NULL ||
tmp_lat1d == NULL || tmp_lon1d == NULL || tmp_fi == NULL) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"rcm2points: Unable to coerce input lat/lon arrays to double precision");
return(NhlFATAL);
}
/*
* Force ncrit to 1, since it is not used yet.
*/
tmp_ncrit = 1;
/*
* Check if any attributes have been attached to opt.
*/
set_search_width = False;
stack_entry = _NclGetArg(5, 6, DONT_CARE);
switch (stack_entry.kind) {;;
case NclStk_VAR:
if (stack_entry.u.data_var->var.att_id != -1) {;;
attr_obj = (NclAtt) _NclGetObj(stack_entry.u.data_var->var.att_id);
if (attr_obj == NULL) {;
break;
};
}
else {
/*
* att_id == -1 ==> no optional args given.
*/
break;
}
/*
* Get optional arguments.
*/
if (attr_obj->att.n_atts > 0) {
/*
* Get list of attributes.
*/
attr_list = attr_obj->att.att_list;
/*
* Loop through attributes and check them. The current ones recognized are:
*
* "search_width"
*/
while (attr_list != NULL) {
if (!strcmp(attr_list->attname, "search_width")) {
if(attr_list->attvalue->multidval.data_type != NCL_int) {
NhlPError(NhlWARNING,NhlEUNKNOWN,"rcm2points: The 'search_width' attribute must be an integer, defaulting to 1.");
search_width = 1;
}
else {
search_width = *(int*) attr_list->attvalue->multidval.val;
if(search_width <= 0) {
NhlPError(NhlWARNING,NhlEUNKNOWN,"rcm2points: The 'search_width' attribute is < 0. Defaulting to 1.");
search_width = 1;
}
else {
set_search_width = True;
}
}
}
attr_list = attr_list->next;
}
}
default:
break;
}
/*
* If user didn't set search_width, then set it here.
*/
if(!set_search_width) search_width = 1;
NGCALLF(drcm2points,DRCM2POINTS)(&ingrid,&inlat2d,&inlon2d,tmp_lat2d,
tmp_lon2d,tmp_fi,&inlat1d,tmp_lat1d,
tmp_lon1d,tmp_fo,&missing_dfi.doubleval,
opt,&tmp_ncrit,&search_width,&ier);
if(ier) {
if(ier == 1) {
NhlPError(NhlWARNING,NhlEUNKNOWN,"rcm2points: not enough points in input/output array");
}
if(2 <= ier && ier <= 5) {
NhlPError(NhlWARNING,NhlEUNKNOWN,"rcm2points: lat2d, lon2d, lat1d, lon1d must be monotonically increasing");
}
set_subset_output_missing(fo,0,type_fo,size_fo,missing_dfi.doubleval);
}
else {
if(type_fo != NCL_double) {
coerce_output_float_only(fo,tmp_fo,size_fo,0);
}
}
/*
* Free temp arrays.
*/
if(type_lat2d != NCL_double) NclFree(tmp_lat2d);
if(type_lon2d != NCL_double) NclFree(tmp_lon2d);
if(type_lat1d != NCL_double) NclFree(tmp_lat1d);
if(type_lon1d != NCL_double) NclFree(tmp_lon1d);
if(type_fi != NCL_double) NclFree(tmp_fi);
if(type_fo != NCL_double) NclFree(tmp_fo);
/*
* Return.
*/
ret = NclReturnValue(fo,ndims_fo,dsizes_fo,&missing_fo,type_fo,0);
NclFree(dsizes_fo);
return(ret);
}

612
ncl_reference/rip_cape.f
Unescape View File

@ -1,612 +0,0 @@
c======================================================================
c
c !IROUTINE: capecalc3d -- Calculate CAPE and CIN
c
c !DESCRIPTION:
c
c If i3dflag=1, this routine calculates CAPE and CIN (in m**2/s**2,
c or J/kg) for every grid point in the entire 3D domain (treating
c each grid point as a parcel). If i3dflag=0, then it
c calculates CAPE and CIN only for the parcel with max theta-e in
c the column, (i.e. something akin to Colman's MCAPE). By "parcel",
c we mean a 500-m deep parcel, with actual temperature and moisture
c averaged over that depth.
c
c In the case of i3dflag=0,
c CAPE and CIN are 2D fields that are placed in the k=mkzh slabs of
c the cape and cin arrays. Also, if i3dflag=0, LCL and LFC heights
c are put in the k=mkzh-1 and k=mkzh-2 slabs of the cin array.
c
c ASSUMPTIONS:
c
c !REVISION HISTORY:
c 2005-May-15 - Mark T. Stoelinga - oringinal version from RIP4
c 2005-Nov-28 - J. Schramm - modified to run outside of RIP4 with
c 2012-Jul-18 - M. Haley - modified to change/add missing value.
c NCL
c
c !INTERFACE:
c ------------------------------------------------------------------
C NCLFORTSTART
SUBROUTINE DCAPECALC3D(PRS,TMK,QVP,GHT,TER,SFP,CAPE,CIN,CMSG,
+ MIY,MJX,MKZH,I3DFLAG,TER_FOLLOW,PSAFILE)
c
IMPLICIT NONE
INTEGER MIY,MJX,MKZH,I3DFLAG,TER_FOLLOW
DOUBLE PRECISION PRS(MIY,MJX,MKZH)
DOUBLE PRECISION TMK(MIY,MJX,MKZH)
DOUBLE PRECISION QVP(MIY,MJX,MKZH)
DOUBLE PRECISION GHT(MIY,MJX,MKZH)
DOUBLE PRECISION TER(MIY,MJX)
DOUBLE PRECISION SFP(MIY,MJX)
DOUBLE PRECISION CAPE(MIY,MJX,MKZH)
DOUBLE PRECISION CIN(MIY,MJX,MKZH)
DOUBLE PRECISION CMSG
CHARACTER*(*) PSAFILE
C NCLEND
c Local variables
INTEGER I,J,K,ILCL,IUP,KEL,KK,KLCL,KLEV,KLFC,KMAX,KPAR,KPAR1,KPAR2
DOUBLE PRECISION DAVG,ETHMAX,Q,T,P,E,ETH,TLCL,ZLCL
DOUBLE PRECISION CP,EPS,GAMMA,GAMMAMD,RGAS,RGASMD,TLCLC1,TLCLC2,
+ TLCLC3,TLCLC4
DOUBLE PRECISION CPMD,THTECON1,THTECON2,THTECON3
DOUBLE PRECISION CELKEL,EZERO,ESLCON1,ESLCON2,GRAV
DOUBLE PRECISION PAVG,VIRTUAL,P1,P2,PP1,PP2,TH,TOTTHE,TOTQVP,
+ TOTPRS
DOUBLE PRECISION CPM,DELTAP,ETHPARI,GAMMAM,GHTPARI,QVPPARI,
+ PRSPARI,TMKPARI
DOUBLE PRECISION FACDEN,FAC1,FAC2,QVPLIFT,TMKLIFT,TVENV,TVLIFT,
+ GHTLIFT
DOUBLE PRECISION ESLIFT,TMKENV,QVPENV,TONPSADIABAT
DOUBLE PRECISION BENAMIN,DZ,PUP,PDN
DOUBLE PRECISION BUOY(150),ZREL(150),BENACCUM(150),
+ PRSF(MIY,MJX,MKZH)
DOUBLE PRECISION PSADITHTE(150),PSADIPRS(150),PSADITMK(150,150)
c
C The comments were taken from a Mark Stoelinga email, 23 Apr 2007,
C in response to a user getting the "Outside of lookup table bounds"
C error message.
C
C TMKPARI - Initial temperature of parcel, K
C Values of 300 okay. (Not sure how much from this you can stray.)
C
C PRSPARI - Initial pressure of parcel, hPa
C Values of 980 okay. (Not sure how much from this you can stray.)
C
C THTECON1, THTECON2, THTECON3
C These are all constants, the first in K and the other two have
C no units. Values of 3376, 2.54, and 0.81 were stated as being
C okay.
C
C TLCL - The temperature at the parcel's lifted condensation level, K
C should be a reasonable atmospheric temperature around 250-300 K
C (398 is "way too high")
C
C QVPPARI - The initial water vapor mixing ratio of the parcel,
C kg/kg (should range from 0.000 to 0.025)
C
c Constants
IUP = 6
CELKEL = 273.15D0
GRAV = 9.81D0
C hPa
EZERO = 6.112D0
ESLCON1 = 17.67D0
ESLCON2 = 29.65D0
EPS = 0.622D0
C J/K/kg
RGAS = 287.04D0
C J/K/kg Note: not using Bolton's value of 1005.7
CP = 1004.D0
GAMMA = RGAS/CP
C cp_moist=cp*(1.+cpmd*qvp)
CPMD = .887D0
C rgas_moist=rgas*(1.+rgasmd*qvp)
RGASMD = .608D0
C gamma_moist=gamma*(1.+gammamd*qvp)
GAMMAMD = RGASMD - CPMD
TLCLC1 = 2840.D0
TLCLC2 = 3.5D0
TLCLC3 = 4.805D0
TLCLC4 = 55.D0
C K
THTECON1 = 3376.D0
THTECON2 = 2.54D0
THTECON3 = .81D0
c
c Calculated the pressure at full sigma levels (a set of pressure
c levels that bound the layers represented by the vertical grid points)
CALL DPFCALC(PRS,SFP,PRSF,MIY,MJX,MKZH,TER_FOLLOW)
c
c Before looping, set lookup table for getting temperature on
c a pseudoadiabat.
c
CALL DLOOKUP_TABLE(PSADITHTE,PSADIPRS,PSADITMK,PSAFILE)
c
C do j=1,mjx-1
DO J = 1,MJX
C do i=1,miy-1
DO I = 1,MIY
CAPE(I,J,1) = 0.D0
CIN(I,J,1) = 0.D0
c
IF (I3DFLAG.EQ.1) THEN
KPAR1 = 2
KPAR2 = MKZH
ELSE
c
c Find parcel with max theta-e in lowest 3 km AGL.
c
ETHMAX = -1.D0
DO K = MKZH,1,-1
IF (GHT(I,J,K)-TER(I,J).LT.3000.D0) THEN
Q = MAX(QVP(I,J,K),1.D-15)
T = TMK(I,J,K)
P = PRS(I,J,K)
E = Q*P/ (EPS+Q)
TLCL = TLCLC1/ (LOG(T**TLCLC2/E)-TLCLC3) +
+ TLCLC4
ETH = T* (1000.D0/P)**
+ (GAMMA* (1.D0+GAMMAMD*Q))*
+ EXP((THTECON1/TLCL-THTECON2)*Q*
+ (1.D0+THTECON3*Q))
IF (ETH.GT.ETHMAX) THEN
KLEV = K
ETHMAX = ETH
END IF
END IF
END DO
KPAR1 = KLEV
KPAR2 = KLEV
c
c Establish average properties of that parcel
c (over depth of approximately davg meters)
c
c davg=.1
DAVG = 500.D0
PAVG = DAVG*PRS(I,J,KPAR1)*GRAV/
+ (RGAS*VIRTUAL(TMK(I,J,KPAR1),QVP(I,J,KPAR1)))
P2 = MIN(PRS(I,J,KPAR1)+.5D0*PAVG,PRSF(I,J,MKZH))
P1 = P2 - PAVG
TOTTHE = 0.D0
TOTQVP = 0.D0
TOTPRS = 0.D0
DO K = MKZH,2,-1
IF (PRSF(I,J,K).LE.P1) GO TO 35
IF (PRSF(I,J,K-1).GE.P2) GO TO 34
P = PRS(I,J,K)
PUP = PRSF(I,J,K)
PDN = PRSF(I,J,K-1)
Q = MAX(QVP(I,J,K),1.D-15)
TH = TMK(I,J,K)* (1000.D0/P)**
+ (GAMMA* (1.D0+GAMMAMD*Q))
PP1 = MAX(P1,PDN)
PP2 = MIN(P2,PUP)
IF (PP2.GT.PP1) THEN
DELTAP = PP2 - PP1
TOTQVP = TOTQVP + Q*DELTAP
TOTTHE = TOTTHE + TH*DELTAP
TOTPRS = TOTPRS + DELTAP
END IF
34 CONTINUE
END DO
35 CONTINUE
QVPPARI = TOTQVP/TOTPRS
TMKPARI = (TOTTHE/TOTPRS)*
+ (PRS(I,J,KPAR1)/1000.D0)** (GAMMA*
+ (1.D0+GAMMAMD*QVP(I,J,KPAR1)))
END IF
c
DO KPAR = KPAR1,KPAR2
c
c Calculate temperature and moisture properties of parcel
c (Note, qvppari and tmkpari already calculated above for 2D case.)
c
IF (I3DFLAG.EQ.1) THEN
QVPPARI = QVP(I,J,KPAR)
TMKPARI = TMK(I,J,KPAR)
END IF
PRSPARI = PRS(I,J,KPAR)
GHTPARI = GHT(I,J,KPAR)
GAMMAM = GAMMA* (1.D0+GAMMAMD*QVPPARI)
CPM = CP* (1.D0+CPMD*QVPPARI)
c
E = MAX(1.D-20,QVPPARI*PRSPARI/ (EPS+QVPPARI))
TLCL = TLCLC1/ (LOG(TMKPARI**TLCLC2/E)-TLCLC3) +
+ TLCLC4
ETHPARI = TMKPARI* (1000.D0/PRSPARI)**
+ (GAMMA* (1.D0+GAMMAMD*QVPPARI))*
+ EXP((THTECON1/TLCL-THTECON2)*QVPPARI*
+ (1.D0+THTECON3*QVPPARI))
ZLCL = GHTPARI + (TMKPARI-TLCL)/ (GRAV/CPM)
c
c Calculate buoyancy and relative height of lifted parcel at
c all levels, and store in bottom up arrays. Add a level at the LCL,
c and at all points where buoyancy is zero.
c
C for arrays that go bottom to top
KK = 0
ILCL = 0
IF (GHTPARI.GE.ZLCL) THEN
c
c initial parcel already saturated or supersaturated.
c
ILCL = 2
KLCL = 1
END IF
DO K = KPAR,1,-1
C for arrays that go bottom to top
33 KK = KK + 1
C model level is below LCL
IF (GHT(I,J,K).LT.ZLCL) THEN
QVPLIFT = QVPPARI
TMKLIFT = TMKPARI - GRAV/CPM*
+ (GHT(I,J,K)-GHTPARI)
TVENV = VIRTUAL(TMK(I,J,K),QVP(I,J,K))
TVLIFT = VIRTUAL(TMKLIFT,QVPLIFT)
GHTLIFT = GHT(I,J,K)
ELSE IF (GHT(I,J,K).GE.ZLCL .AND. ILCL.EQ.0) THEN
c
c This model level and previous model level straddle the LCL,
c so first create a new level in the bottom-up array, at the LCL.
c
TMKLIFT = TLCL
QVPLIFT = QVPPARI
FACDEN = GHT(I,J,K) - GHT(I,J,K+1)
FAC1 = (ZLCL-GHT(I,J,K+1))/FACDEN
FAC2 = (GHT(I,J,K)-ZLCL)/FACDEN
TMKENV = TMK(I,J,K+1)*FAC2 + TMK(I,J,K)*FAC1
QVPENV = QVP(I,J,K+1)*FAC2 + QVP(I,J,K)*FAC1
TVENV = VIRTUAL(TMKENV,QVPENV)
TVLIFT = VIRTUAL(TMKLIFT,QVPLIFT)
GHTLIFT = ZLCL
ILCL = 1
ELSE
TMKLIFT = TONPSADIABAT(ETHPARI,PRS(I,J,K),
+ PSADITHTE,PSADIPRS,PSADITMK,GAMMA)
ESLIFT = EZERO*EXP(ESLCON1* (TMKLIFT-CELKEL)/
+ (TMKLIFT-ESLCON2))
QVPLIFT = EPS*ESLIFT/ (PRS(I,J,K)-ESLIFT)
TVENV = VIRTUAL(TMK(I,J,K),QVP(I,J,K))
TVLIFT = VIRTUAL(TMKLIFT,QVPLIFT)
GHTLIFT = GHT(I,J,K)
END IF
C buoyancy
BUOY(KK) = GRAV* (TVLIFT-TVENV)/TVENV
ZREL(KK) = GHTLIFT - GHTPARI
IF ((KK.GT.1).AND.
+ (BUOY(KK)*BUOY(KK-1).LT.0.0D0)) THEN
c
c Parcel ascent curve crosses sounding curve, so create a new level
c in the bottom-up array at the crossing.
c
KK = KK + 1
BUOY(KK) = BUOY(KK-1)
ZREL(KK) = ZREL(KK-1)
BUOY(KK-1) = 0.D0
ZREL(KK-1) = ZREL(KK-2) +
+ BUOY(KK-2)/ (BUOY(KK-2)-
+ BUOY(KK))* (ZREL(KK)-ZREL(KK-2))
END IF
IF (ILCL.EQ.1) THEN
KLCL = KK
ILCL = 2
GO TO 33
END IF
END DO
KMAX = KK
IF (KMAX.GT.150) THEN
print *,
+ 'capecalc3d: kmax got too big. kmax=',KMAX
STOP
END IF
c
c If no LCL was found, set klcl to kmax. It is probably not really
c at kmax, but this will make the rest of the routine behave
c properly.
c
IF (ILCL.EQ.0) KLCL=KMAX
c
c Get the accumulated buoyant energy from the parcel's starting
c point, at all levels up to the top level.
c
BENACCUM(1) = 0.0D0
BENAMIN = 9D9
DO K = 2,KMAX
DZ = ZREL(K) - ZREL(K-1)
BENACCUM(K) = BENACCUM(K-1) +
+ .5D0*DZ* (BUOY(K-1)+BUOY(K))
IF (BENACCUM(K).LT.BENAMIN) THEN
BENAMIN = BENACCUM(K)
END IF
END DO
c
c Determine equilibrium level (EL), which we define as the highest
c level of non-negative buoyancy above the LCL. Note, this may be
c the top level if the parcel is still buoyant there.
c
DO K = KMAX,KLCL,-1
IF (BUOY(K).GE.0.D0) THEN
C k of equilibrium level
KEL = K
GO TO 50
END IF
END DO
c
c If we got through that loop, then there is no non-negative
c buoyancy above the LCL in the sounding. In these situations,
c both CAPE and CIN will be set to -0.1 J/kg. (See below about
c missing values in V6.1.0). Also, where CAPE is
c non-zero, CAPE and CIN will be set to a minimum of +0.1 J/kg, so
c that the zero contour in either the CIN or CAPE fields will
c circumscribe regions of non-zero CAPE.
c
c In V6.1.0 of NCL, we added a _FillValue attribute to the return
c value of this function. At that time we decided to change -0.1
c to a more appropriate missing value, which is passed into this
c routine as CMSG.
c
c CAPE(I,J,KPAR) = -0.1D0
c CIN(I,J,KPAR) = -0.1D0
CAPE(I,J,KPAR) = CMSG
CIN(I,J,KPAR) = CMSG
KLFC = KMAX
c
GO TO 102
c
50 CONTINUE
c
c If there is an equilibrium level, then CAPE is positive. We'll
c define the level of free convection (LFC) as the point below the
c EL, but at or above the LCL, where accumulated buoyant energy is a
c minimum. The net positive area (accumulated buoyant energy) from
c the LFC up to the EL will be defined as the CAPE, and the net
c negative area (negative of accumulated buoyant energy) from the
c parcel starting point to the LFC will be defined as the convective
c inhibition (CIN).
c
c First get the LFC according to the above definition.
c
BENAMIN = 9D9
KLFC = KMAX
DO K = KLCL,KEL
IF (BENACCUM(K).LT.BENAMIN) THEN
BENAMIN = BENACCUM(K)
KLFC = K
END IF
END DO
c
c Now we can assign values to cape and cin
c
CAPE(I,J,KPAR) = MAX(BENACCUM(KEL)-BENAMIN,0.1D0)
CIN(I,J,KPAR) = MAX(-BENAMIN,0.1D0)
c
c CIN is uninteresting when CAPE is small (< 100 J/kg), so set
c CIN to -0.1 (see note about missing values in V6.1.0) in
c that case.
c
c In V6.1.0 of NCL, we added a _FillValue attribute to the return
c value of this function. At that time we decided to change -0.1
c to a more appropriate missing value, which is passed into this
c routine as CMSG.
c
C IF (CAPE(I,J,KPAR).LT.100.D0) CIN(I,J,KPAR) = -0.1D0
IF (CAPE(I,J,KPAR).LT.100.D0) CIN(I,J,KPAR) = CMSG
102 CONTINUE
c
END DO
c
IF (I3DFLAG.EQ.0) THEN
CAPE(I,J,MKZH) = CAPE(I,J,KPAR1)
CIN(I,J,MKZH) = CIN(I,J,KPAR1)
C meters AGL
CIN(I,J,MKZH-1) = ZREL(KLCL) + GHTPARI - TER(I,J)
C meters AGL
CIN(I,J,MKZH-2) = ZREL(KLFC) + GHTPARI - TER(I,J)
END IF
c
END DO
END DO
c
RETURN
END
c c
c*********************************************************************c
c c
C NCLFORTSTART
DOUBLE PRECISION FUNCTION TONPSADIABAT(THTE,PRS,PSADITHTE,
& PSADIPRS,PSADITMK,GAMMA)
IMPLICIT NONE
DOUBLE PRECISION THTE
DOUBLE PRECISION PRS
DOUBLE PRECISION PSADITHTE
DOUBLE PRECISION PSADIPRS
DOUBLE PRECISION PSADITMK
DOUBLE PRECISION GAMMA
C NCLEND
DOUBLE PRECISION FRACJT
DOUBLE PRECISION FRACJT2
DOUBLE PRECISION FRACIP
DOUBLE PRECISION FRACIP2
DIMENSION PSADITHTE(150),PSADIPRS(150),PSADITMK(150,150)
INTEGER IP, IPCH, JT, JTCH
c c
c This function gives the temperature (in K) on a moist adiabat
c (specified by thte in K) given pressure in hPa. It uses a
c lookup table, with data that was generated by the Bolton (1980)
c formula for theta_e.
c
c First check if pressure is less than min pressure in lookup table.
c If it is, assume parcel is so dry that the given theta-e value can
c be interpretted as theta, and get temperature from the simple dry
c theta formula.
c
IF (PRS.LE.PSADIPRS(150)) THEN
TONPSADIABAT = THTE* (PRS/1000.D0)**GAMMA
RETURN
END IF
c
c Otherwise, look for the given thte/prs point in the lookup table.
c
DO JTCH = 1,150 - 1
IF (THTE.GE.PSADITHTE(JTCH) .AND.
+ THTE.LT.PSADITHTE(JTCH+1)) THEN
JT = JTCH
GO TO 213
END IF
END DO
JT = -1
213 CONTINUE
DO IPCH = 1,150 - 1
IF (PRS.LE.PSADIPRS(IPCH) .AND. PRS.GT.PSADIPRS(IPCH+1)) THEN
IP = IPCH
GO TO 215
END IF
END DO
IP = -1
215 CONTINUE
IF (JT.EQ.-1 .OR. IP.EQ.-1) THEN
print *,'capecalc3d: ',
+ 'Outside of lookup table bounds. prs,thte=',
+ PRS,THTE
STOP
END IF
FRACJT = (THTE-PSADITHTE(JT))/ (PSADITHTE(JT+1)-PSADITHTE(JT))
FRACJT2 = 1.D0 - FRACJT
FRACIP = (PSADIPRS(IP)-PRS)/ (PSADIPRS(IP)-PSADIPRS(IP+1))
FRACIP2 = 1.D0 - FRACIP
IF (PSADITMK(IP,JT).GT.1D9 .OR. PSADITMK(IP+1,JT).GT.1D9 .OR.
+ PSADITMK(IP,JT+1).GT.1D9 .OR. PSADITMK(IP+1,JT+1).GT.1D9) THEN
print *,'capecalc3d: ',
+ 'Tried to access missing temperature in lookup table.',
+ 'Prs and Thte probably unreasonable. prs,thte=',PRS,THTE
STOP
END IF
TONPSADIABAT = FRACIP2*FRACJT2*PSADITMK(IP,JT) +
+ FRACIP*FRACJT2*PSADITMK(IP+1,JT) +
+ FRACIP2*FRACJT*PSADITMK(IP,JT+1) +
+ FRACIP*FRACJT*PSADITMK(IP+1,JT+1)
c
RETURN
END
c c
c*********************************************************************c
SUBROUTINE DLOOKUP_TABLE(PSADITHTE,PSADIPRS,PSADITMK,FNAME)
DOUBLE PRECISION PSADITHTE
DOUBLE PRECISION PSADIPRS
DOUBLE PRECISION PSADITMK
c Set up lookup table for getting temperature on a pseudoadiabat.
c (Borrow the unit number for the stationlist, just for the moment.)
c
C CHARACTER*15 FNAME
CHARACTER*(*) FNAME
DIMENSION PSADITHTE(150),PSADIPRS(150),PSADITMK(150,150)
C FNAME = 'psadilookup.dat'
IUSTNLIST = 33
OPEN (UNIT=IUSTNLIST,FILE=FNAME,FORM='formatted',STATUS='old')
DO I = 1,14
READ (IUSTNLIST,FMT=*)
END DO
READ (IUSTNLIST,FMT=*) NTHTE,NPRS
IF (NTHTE.NE.150 .OR. NPRS.NE.150) THEN
WRITE (IUP,FMT=*)
+ 'Number of pressure or theta_e levels in lookup table'
WRITE (IUP,FMT=*) 'file not = 150. Check lookup table file.'
STOP
END IF
READ (IUSTNLIST,FMT=173) (PSADITHTE(JT),JT=1,NTHTE)
READ (IUSTNLIST,FMT=173) (PSADIPRS(IP),IP=1,NPRS)
READ (IUSTNLIST,FMT=173) ((PSADITMK(IP,JT),IP=1,NPRS),JT=1,NTHTE)
173 FORMAT (5D15.7)
CLOSE (IUSTNLIST)
RETURN
END
c c
c*********************************************************************c
c c
SUBROUTINE DPFCALC(PRS,SFP,PF,MIY,MJX,MKZH,TER_FOLLOW)
DOUBLE PRECISION PRS
DOUBLE PRECISION SFP
DOUBLE PRECISION PF
c
c Historically, this routine calculated the pressure at full sigma
c levels when RIP was specifically designed for MM4/MM5 output.
c With the new generalized RIP (Feb '02), this routine is still
c intended to calculate a set of pressure levels that bound the
c layers represented by the vertical grid points, although no such
c layer boundaries are assumed to be defined. The routine simply
c uses the midpoint between the pressures of the vertical grid
c points as the bounding levels. The array only contains mkzh
c levels, so the pressure of the top of the uppermost layer is
c actually excluded. The kth value of pf is the lower bounding
c pressure for the layer represented by kth data level. At the
c lower bounding level of the lowest model layer, it uses the
c surface pressure, unless the data set is pressure-level data, in
c which case it assumes the lower bounding pressure level is as far
c below the lowest vertical level as the upper bounding pressure
c level is above.
c
DIMENSION PRS(MIY,MJX,MKZH),SFP(MIY,MJX),PF(MIY,MJX,MKZH)
INTEGER TER_FOLLOW
c
C do j=1,mjx-1 Artifact of MM5
DO J = 1,MJX
C do i=1,miy-1 staggered grid
DO I = 1,MIY
DO K = 1,MKZH
IF (K.EQ.MKZH) THEN
C terrain-following data
IF (TER_FOLLOW.EQ.1) THEN
PF(I,J,K) = SFP(I,J)
C pressure-level data
ELSE
PF(I,J,K) = .5D0* (3.D0*PRS(I,J,K)-
+ PRS(I,J,K-1))
END IF
ELSE
PF(I,J,K) = .5D0* (PRS(I,J,K+1)+PRS(I,J,K))
END IF
END DO
END DO
END DO
c
RETURN
END
c======================================================================
c
c !IROUTINE: VIRTUAL -- Calculate virtual temperature (K)
c
c !DESCRIPTION:
c
c This function returns a single value of virtual temperature in
c K, given temperature in K and mixing ratio in kg/kg. For an
c array of virtual temperatures, use subroutine VIRTUAL_TEMP.
c
c !INPUT:
c RATMIX - water vapor mixing ratio (kg/kg)
c TEMP - temperature (K)
c
c !OUTPUT:
c TV - Virtual temperature (K)
c
c !ASSUMPTIONS:
c
c !REVISION HISTORY:
c 2009-March - Mark T. Stoelinga - from RIP4.5
c 2010-August - J. Schramm - modified to run with NCL and ARW wrf output
c
c ------------------------------------------------------------------
C NCLFORTSTART
DOUBLE PRECISION FUNCTION VIRTUAL(TEMP,RATMIX)
IMPLICIT NONE
DOUBLE PRECISION TEMP,RATMIX
C NCLEND
DOUBLE PRECISION EPS
EPS = 0.622D0
VIRTUAL = TEMP* (EPS+RATMIX)/ (EPS* (1.D0+RATMIX))
RETURN
END

13253
ncl_reference/wrfW.c
View File

File diff suppressed because it is too large Load Diff

402
ncl_reference/wrf_bint3d.f
Unescape View File

@ -1,402 +0,0 @@
C
C premaptform.f and maptform.f copied from RIP/src
C By So-Young Ha on Sep 29, 2005.
C
C
C NCLFORTSTART
SUBROUTINE DMAPTFORM(DSKMC,MIYCORS,MJXCORS,NPROJ,XLATC,XLONC,
+ TRUE1,TRUE2,RIY,RJX,RLAT,RLON,IDIR)
C
C Input vars: DSKMC, MIYCORS, MJXCORS, NPROJ, XLATC, XLONC,
C NPROJ, IDIR
C Input/output vars: RIY, RIX, RLAT
C Output vars: TRUE1, TRUE2, RLON
C
C
C Possible NCL interface:
C
C wrf_maptform(dskmc, miycors, mjxcors, nproj, xlatc, xlonc, riy, rjx,
C idir, rlat, rlon, opts)
C
C where opts could contain the TRUE1 and TRUE2 information in some fashion.
C
DOUBLE PRECISION PI_MPTF
DOUBLE PRECISION RPD_MPTF
DOUBLE PRECISION REARTH_MPTF
DOUBLE PRECISION DSKMC_MPTF
DOUBLE PRECISION XLONC_MPTF
DOUBLE PRECISION CIY_MPTF
DOUBLE PRECISION CJX_MPTF
DOUBLE PRECISION CONE_MPTF
DOUBLE PRECISION CONEI_MPTF
DOUBLE PRECISION C1_MPTF
DOUBLE PRECISION C2_MPTF
DOUBLE PRECISION YC_MPTF
DOUBLE PRECISION COTRUE1
DOUBLE PRECISION YPOINT
DOUBLE PRECISION XPOINT
DOUBLE PRECISION DLON
C
c This routine converts a coarse domain dot grid point, <riy,rjx>,
c into a lat/lon point <rlat,rlon> if idir=1, or vice versa if
c idir=-1. It works for Lambert Conformal (LC,1),
c Polar Stereographic (ST,2), or Mercator (ME,3) projections,
c with any true latitide(s).
c It is assumed that premaptform has been called prior to this so
c that the proper constants have been placed in the common block
c called mptf, which should be declared in (and only in) the
c main program and routines maptform (this routine) and premaptform.
c
C Input, Output Args
INTEGER MIYCORS,MJXCORS,NPROJ
DOUBLE PRECISION DSKMC,XLATC,XLONC,TRUE1,TRUE2
INTEGER IDIR
C Latitude (-90->90 deg N)
DOUBLE PRECISION RLAT
C Longitude (-180->180 E)
DOUBLE PRECISION RLON
C Cartesian X coordinate
DOUBLE PRECISION RIY
C Cartesian Y coordinate
DOUBLE PRECISION RJX
C NCLEND
c ===========
c premaptform
c ===========
C 3.1415...
PI_MPTF = 4.D0*ATAN(1.D0)
C radians per degree
RPD_MPTF = PI_MPTF/180.D0
C radius of planet, in km
REARTH_MPTF = 6370.949D0
DSKMC_MPTF = DSKMC
XLONC_MPTF = XLONC
NPROJ_MPTF = NPROJ
CIY_MPTF = .5D0* (1.D0+MIYCORS)
CJX_MPTF = .5D0* (1.D0+MJXCORS)
c
C Mercator
IF (NPROJ_MPTF.EQ.3) THEN
c
TRUE1 = 0.D0
TRUE2 = 0.D0
IHM_MPTF = 1
CONE_MPTF = 1.D0
CONEI_MPTF = 1.D0
C1_MPTF = 1.D0
C2_MPTF = 1.D0
YC_MPTF = REARTH_MPTF*LOG((1.D0+SIN(RPD_MPTF*XLATC))/
+ COS(RPD_MPTF*XLATC))
c
C Lambert Comformal or Polar Stereographic
ELSE
c
c Make sure xlatc, true1, and true2 are all in same hemisphere,
c and calculate ihm_mptf.
c
IF (XLATC.GT.0.D0 .AND. TRUE1.GT.0.D0 .AND.
+ TRUE2.GT.0.D0) THEN
IHM_MPTF = 1
ELSE IF (XLATC.LT.0.D0 .AND. TRUE1.LT.0.D0 .AND.
+ TRUE2.LT.0.D0) THEN
IHM_MPTF = -1
ELSE
WRITE (*,FMT=*) 'Invalid latitude parameters for map.'
STOP
END IF
c
c Calculate cone factor
c
IF (NPROJ_MPTF.EQ.1) THEN
IF (TRUE1.NE.TRUE2) THEN
CONE_MPTF = LOG10(COS(RPD_MPTF*TRUE1)/
+ COS(RPD_MPTF*TRUE2))/
+ LOG10(TAN(.25D0*PI_MPTF-
+ IHM_MPTF*.5D0*RPD_MPTF*TRUE1)/
+ TAN(.25D0*PI_MPTF-IHM_MPTF*.5D0*RPD_MPTF*
+ TRUE2))
ELSE
CONE_MPTF = COS(RPD_MPTF* (90.D0-IHM_MPTF*TRUE1))
END IF
ELSE IF (NPROJ_MPTF.EQ.2) THEN
CONE_MPTF = 1.D0
END IF
c
c Calculate other constants
c
CONEI_MPTF = 1.D0/CONE_MPTF
COTRUE1 = IHM_MPTF*90.D0 - TRUE1
IF (NPROJ_MPTF.EQ.1) THEN
C1_MPTF = REARTH_MPTF*SIN(RPD_MPTF*COTRUE1)/
+ (CONE_MPTF* (IHM_MPTF*TAN(.5D0*RPD_MPTF*
+ COTRUE1))**CONE_MPTF)
C2_MPTF = TAN(.5D0*RPD_MPTF*COTRUE1)*
+ (CONE_MPTF/ (IHM_MPTF*REARTH_MPTF*SIN(RPD_MPTF*
+ COTRUE1)))**CONEI_MPTF
YC_MPTF = -C1_MPTF* (IHM_MPTF*
+ TAN(.25D0* (IHM_MPTF*PI_MPTF-
+ 2.D0*RPD_MPTF*XLATC)))**CONE_MPTF
ELSE IF (NPROJ_MPTF.EQ.2) THEN
C1_MPTF = 1.D0 + COS(RPD_MPTF*COTRUE1)
C2_MPTF = 1.D0
YC_MPTF = -REARTH_MPTF*SIN(.5D0*IHM_MPTF*PI_MPTF-
+ RPD_MPTF*XLATC)*C1_MPTF/
+ (1.D0+COS(.5D0*IHM_MPTF*PI_MPTF-RPD_MPTF*XLATC))
END IF
c
END IF
c ========
c maptform
c ========
IF (RLAT.EQ.-90.D0) PRINT *,'maptform:',RIY,RJX,RLAT,RLON,IDIR
C First, deal with idir=1
IF (IDIR.EQ.1) THEN
c
YPOINT = (RIY-CIY_MPTF)*DSKMC_MPTF + YC_MPTF
XPOINT = (RJX-CJX_MPTF)*DSKMC_MPTF
c
IF (NPROJ_MPTF.EQ.3) THEN
RLAT = (2.D0*ATAN(EXP(YPOINT/REARTH_MPTF))-.5D0*PI_MPTF)/
+ RPD_MPTF
RLON = XLONC_MPTF + (XPOINT/REARTH_MPTF)/RPD_MPTF
ELSE IF (NPROJ_MPTF.EQ.1) THEN
RLAT = (.5D0*IHM_MPTF*PI_MPTF-
+ 2.D0*ATAN(C2_MPTF* (SQRT(XPOINT**2+
+ YPOINT**2))**CONEI_MPTF))/RPD_MPTF
RLON = XLONC_MPTF + (CONEI_MPTF*
+ ATAN2(XPOINT,-IHM_MPTF*YPOINT))/RPD_MPTF
ELSE IF (NPROJ_MPTF.EQ.2) THEN
RLAT = (.5D0*IHM_MPTF*PI_MPTF-
+ IHM_MPTF*2.D0*ATAN(SQRT(XPOINT**2+
+ YPOINT**2)/ (REARTH_MPTF*C1_MPTF)))/RPD_MPTF
IF (XPOINT.EQ.0.D0 .AND. YPOINT.EQ.0.D0) THEN
RLON = XLONC_MPTF
ELSE
RLON = XLONC_MPTF + (ATAN2(XPOINT,-IHM_MPTF*YPOINT))/
+ RPD_MPTF
END IF
END IF
RLON = MOD(RLON+900.D0,360.D0) - 180.D0
c
C Otherwise, deal with idir=-1
ELSE
c
DLON = RLON - XLONC_MPTF
IF (DLON.LT.-180.D0) DLON = DLON + 360
IF (DLON.GT.180.D0) DLON = DLON - 360
IF (NPROJ_MPTF.EQ.3) THEN
YPOINT = REARTH_MPTF*LOG((1.D0+SIN(RPD_MPTF*RLAT))/
+ COS(RPD_MPTF*RLAT))
XPOINT = DLON*RPD_MPTF*REARTH_MPTF
ELSE IF (NPROJ_MPTF.EQ.1) THEN
YPOINT = -C1_MPTF* (IHM_MPTF*
+ TAN(.25D0* (IHM_MPTF*PI_MPTF-2.D0*RPD_MPTF*
+ RLAT)))**CONE_MPTF*COS(CONE_MPTF*RPD_MPTF*DLON)
XPOINT = IHM_MPTF*C1_MPTF* (IHM_MPTF*
+ TAN(.25D0* (IHM_MPTF*PI_MPTF-
+ 2.D0*RPD_MPTF*RLAT)))**CONE_MPTF*
+ SIN(CONE_MPTF*RPD_MPTF*DLON)
ELSE IF (NPROJ_MPTF.EQ.2) THEN
YPOINT = -REARTH_MPTF*SIN(.5D0*IHM_MPTF*PI_MPTF-
+ RPD_MPTF*RLAT)*C1_MPTF/ (1.D0+
+ COS(.5D0*IHM_MPTF*PI_MPTF-RPD_MPTF*RLAT))*
+ COS(RPD_MPTF*DLON)
XPOINT = IHM_MPTF*REARTH_MPTF*
+ SIN(.5D0*IHM_MPTF*PI_MPTF-RPD_MPTF*RLAT)*C1_MPTF/
+ (1.D0+COS(.5D0*IHM_MPTF*PI_MPTF-RPD_MPTF*RLAT))*
+ SIN(RPD_MPTF*DLON)
END IF
RIY = (YPOINT-YC_MPTF)/DSKMC_MPTF + CIY_MPTF
RJX = XPOINT/DSKMC_MPTF + CJX_MPTF
c
END IF
RETURN
END
C********************************************************
C NCLFORTSTART
SUBROUTINE DBINT3D(DATA_OUT,OBSII,OBSJJ,DATA_IN,NX,NY,NZ,NOBSICRS,
+ NOBSJCRS,ICRS,JCRS)
C
C Possible NCL interface:
C
C data_out = wrf_bint3d(data_in,obsii,obsjj,icrs,jcrs)
C
C !!! 1_based_array (cols x rows) in fortran <=> 0_based_array
C (rows x cols) in NCL !!!
C !!! Include K-index to make a 3-D array !!!
C
C INPUT VARIABLES
C ---------------
INTEGER ICRS,JCRS,NX,NY,NZ
INTEGER NOBSJCRS,NOBSICRS
DOUBLE PRECISION OBSII(NOBSICRS,NOBSJCRS)
DOUBLE PRECISION OBSJJ(NOBSICRS,NOBSJCRS)
DOUBLE PRECISION DATA_IN(NX,NY,NZ)
C OUTPUT
C ---------------
DOUBLE PRECISION DATA_OUT(NOBSICRS,NOBSJCRS,NZ)
C NCLEND
C LOCAL
DOUBLE PRECISION OBSI,OBSJ
DOUBLE PRECISION DATA_OBS
C
DO K = 1,NZ
DO J = 1,NOBSJCRS
DO I = 1,NOBSICRS
C grid index in lon
OBSI = OBSII(I,J)
C grid index in lat
OBSJ = OBSJJ(I,J)
DATA_OBS = 0.0D0
CALL DBINT(DATA_OBS,OBSI,OBSJ,DATA_IN(1,1,K),NX,NY,
+ ICRS,JCRS)
DATA_OUT(I,J,K) = DATA_OBS
END DO
END DO
END DO
RETURN
END
SUBROUTINE DBINT(PP,XX,YY,LIST,III,JJJ,ICRS,JCRS)
DOUBLE PRECISION PP
DOUBLE PRECISION X
DOUBLE PRECISION Y
DOUBLE PRECISION A
DOUBLE PRECISION B
DOUBLE PRECISION C
DOUBLE PRECISION D
DOUBLE PRECISION E
DOUBLE PRECISION F
DOUBLE PRECISION G
DOUBLE PRECISION H
DOUBLE PRECISION QQ
C
C --- BI-LINEAR INTERPOLATION AMONG FOUR GRID VALUES
C
C INPUT : LIST, XX, YY
C OUTPUT: PP
C
INTEGER ICRS,JCRS,III,JJJ
DOUBLE PRECISION XX,YY
DOUBLE PRECISION LIST(III,JJJ),STL(4,4)
C MASS GRID IN WRF (I-> west-east, J-> south-north)
C
IB = III - ICRS
JB = JJJ - JCRS
PP = 0.0D0
N = 0
I = INT(XX+0.00001D0)
J = INT(YY+0.00001D0)
X = XX - I
Y = YY - J
IF ((ABS(X).GT.0.00001D0) .OR. (ABS(Y).GT.0.00001D0)) THEN
C
DO 2 K = 1,4
KK = I + K
DO 2 L = 1,4
STL(K,L) = 0.D0
LL = J + L
IF ((KK.GE.1) .AND. (KK.LE.IB) .AND. (LL.LE.JB) .AND.
+ (LL.GE.1)) THEN
STL(K,L) = LIST(KK,LL)
N = N + 1
C .. a zero value inside the domain being set to 1.E-12:
IF (STL(K,L).EQ.0.D0) STL(K,L) = 1.D-12
END IF
2 CONTINUE
C
CALL DONED(A,X,STL(1,1),STL(2,1),STL(3,1),STL(4,1))
CALL DONED(B,X,STL(1,2),STL(2,2),STL(3,2),STL(4,2))
CALL DONED(C,X,STL(1,3),STL(2,3),STL(3,3),STL(4,3))
CALL DONED(D,X,STL(1,4),STL(2,4),STL(3,4),STL(4,4))
C
C .. CHECK TANGENT LINEAR OF ONED, SAVE BASIC STATE:
C WRITE(20) XX,YY,Y,A,B,C,D
C
CALL DONED(PP,Y,A,B,C,D)
IF (N.NE.16) THEN
CALL DONED(E,Y,STL(1,1),STL(1,2),STL(1,3),STL(1,4))
CALL DONED(F,Y,STL(2,1),STL(2,2),STL(2,3),STL(2,4))
CALL DONED(G,Y,STL(3,1),STL(3,2),STL(3,3),STL(3,4))
CALL DONED(H,Y,STL(4,1),STL(4,2),STL(4,3),STL(4,4))
C .. CHECK TANGENT LINEAR OF ONED, SAVE BASIC STATE:
C WRITE(20) XX,YY,X,E,F,G,H
C
CALL DONED(QQ,X,E,F,G,H)
PP = (PP+QQ)*0.5D0
END IF
C
ELSE
C
PP = LIST(I,J)
END IF
C
RETURN
END
SUBROUTINE DONED(Y,X,A,B,C,D)
DOUBLE PRECISION Y
DOUBLE PRECISION X
DOUBLE PRECISION A
DOUBLE PRECISION B
DOUBLE PRECISION C
DOUBLE PRECISION D
DOUBLE PRECISION ONE
C
C .. Input : X, A, B, C, D
C Output: Y
C 1, 2, 3, and 4 points interpolation:
C In this subroutine, the zero value of A, B, C, D means that
C point outside the domain.
C
C .. 1-point:
C .. take the value at the second point:
IF (X.EQ.0.D0) THEN
ONE = B
C .. take the value at the third point:
ELSE IF (X.EQ.1.D0) THEN
ONE = C
C .. the point X outside the range:
ELSE IF (B*C.EQ.0.D0) THEN
ONE = 0.D0
ELSE
IF (A*D.EQ.0.D0) THEN
C .. 3-point interpolation:
IF (A.NE.0.D0) THEN
ONE = B + X* (0.5D0* (C-A)+X* (0.5D0* (C+A)-B))
ELSE IF (D.NE.0.D0) THEN
ONE = C + (1.0D0-X)* (0.5D0* (B-D)+
+ (1.0D0-X)* (0.5D0* (B+D)-C))
ELSE
C .. 2-point interpolation:
ONE = B* (1.0D0-X) + C*X
END IF
ELSE
C .. 4-point interpolation:
ONE = (1.0D0-X)* (B+X* (0.5D0* (C-A)+X* (0.5D0* (C+A)-B)))
+ + X* (C+ (1.0D0-X)* (0.5D0* (B-D)+ (1.0D0-
+ X)* (0.5D0* (B+D)-C)))
END IF
END IF
C
Y = ONE
C
RETURN
END

763
ncl_reference/wrf_cloud_topW.c
Unescape View File

@ -1,763 +0,0 @@
#include <stdio.h>
#include "wrapper.h"
extern void NGCALLF(wrfcttcalc,WRFCTTCALC)(double *, double *, double *,
double *, double *, double *,
double *, double *, int *,
int *, int *, int *);
extern NclDimRec *get_wrf_dim_info(int,int,int,ng_size_t*);
NhlErrorTypes wrf_ctt_W( void )
{
/*
* Input variables
*/
/*
* Argument # 0
*/
void *pres;
double *tmp_pres;
int ndims_pres;
ng_size_t dsizes_pres[NCL_MAX_DIMENSIONS];
NclBasicDataTypes type_pres;
/*
* Argument # 1
*/
void *tk;
double *tmp_tk;
int ndims_tk;
ng_size_t dsizes_tk[NCL_MAX_DIMENSIONS];
NclBasicDataTypes type_tk;
/*
* Argument # 2
*/
void *qci;
double *tmp_qci;
int ndims_qci;
ng_size_t dsizes_qci[NCL_MAX_DIMENSIONS];
NclBasicDataTypes type_qci;
/*
* Argument # 3
*/
void *qcw;
double *tmp_qcw;
int ndims_qcw;
ng_size_t dsizes_qcw[NCL_MAX_DIMENSIONS];
NclBasicDataTypes type_qcw;
/*
* Argument # 4
*/
void *qvp;
double *tmp_qvp;
int ndims_qvp;
ng_size_t dsizes_qvp[NCL_MAX_DIMENSIONS];
NclBasicDataTypes type_qvp;
/*
* Argument # 5
*/
void *ght;
double *tmp_ght;
int ndims_ght;
ng_size_t dsizes_ght[NCL_MAX_DIMENSIONS];
NclBasicDataTypes type_ght;
/*
* Argument # 6
*/
void *ter;
double *tmp_ter;
int ndims_ter;
ng_size_t dsizes_ter[NCL_MAX_DIMENSIONS];
NclBasicDataTypes type_ter;
/*
* Arguments # 7
*/
int *haveqci;
/*
* Variable for getting/setting dimension name info.
*/
NclDimRec *dim_info = NULL;
NclDimRec *dim_info_ght = NULL;
/*
* Return variable and attributes
*/
void *ctt;
NclQuark *description, *units;
char *cdescription, *cunits;
double *tmp_ctt;
int ndims_ctt;
ng_size_t *dsizes_ctt;
NclBasicDataTypes type_ctt;
NclObjClass type_obj_ctt;
/*
* Various
*/
ng_size_t nlev, nlat, nlon, nlevlatlon, nlatlon;
ng_size_t index_pres, index_ter, index_ctt;
ng_size_t i, size_leftmost, size_output;
int inlev, inlat, inlon;
/*
* Variables for returning the output array with attributes attached.
*/
int att_id;
ng_size_t dsizes[1];
NclMultiDValData att_md, return_md;
NclVar tmp_var;
NclStackEntry return_data;
/*
* Retrieve parameters.
*
* Note any of the pointer parameters can be set to NULL, which
* implies you don't care about its value.
*/
/*
* Get argument # 0
*/
/*
* Get argument # 1
*/
pres = (void*)NclGetArgValue(
0,
8,
&ndims_pres,
dsizes_pres,
NULL,
NULL,
&type_pres,
DONT_CARE);
if(ndims_pres < 3 || ndims_pres > 4) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: The pres array must be 3D or 4D");
return(NhlFATAL);
}
nlev = dsizes_pres[ndims_pres-3];
nlat = dsizes_pres[ndims_pres-2];
nlon = dsizes_pres[ndims_pres-1];
/*
* Test dimension sizes.
*/
if(nlev > INT_MAX || nlat > INT_MAX || nlon > INT_MAX) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: one of bottom_top, south_north, or west_east is greater than INT_MAX");
return(NhlFATAL);
}
inlev = (int) nlev;
inlat = (int) nlat;
inlon = (int) nlon;
/*
* Get argument # 1
*/
tk = (void*)NclGetArgValue(
1,
8,
&ndims_tk,
dsizes_tk,
NULL,
NULL,
&type_tk,
DONT_CARE);
/*
* Check dimension sizes.
*/
if(ndims_tk != ndims_pres) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: The tk and pres arrays must have the same dimensionality");
return(NhlFATAL);
}
else {
for(i = 0; i < ndims_pres; i++) {
if(dsizes_tk[i] != dsizes_pres[i]) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: The tk and pres arrays must have the same dimensionality");
return(NhlFATAL);
}
}
}
/*
* Get argument # 2
*/
qci = (void*)NclGetArgValue(
2,
8,
&ndims_qci,
dsizes_qci,
NULL,
NULL,
&type_qci,
DONT_CARE);
/*
* Check dimension sizes.
*/
if(ndims_qci != ndims_pres) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: The qci and pres arrays must have the same dimensionality");
return(NhlFATAL);
}
else {
for(i = 0; i < ndims_pres; i++) {
if(dsizes_qci[i] != dsizes_pres[i]) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: The qci and pres arrays must have the same dimensionality");
return(NhlFATAL);
}
}
}
/*
* Get argument # 3
*/
qcw = (void*)NclGetArgValue(
3,
8,
&ndims_qcw,
dsizes_qcw,
NULL,
NULL,
&type_qcw,
DONT_CARE);
/*
* Check dimension sizes.
*/
if(ndims_qcw != ndims_pres) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: The qcw and pres arrays must have the same dimensionality");
return(NhlFATAL);
}
else {
for(i = 0; i < ndims_pres; i++) {
if(dsizes_qcw[i] != dsizes_pres[i]) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: The qcw and pres arrays must have the same dimensionality");
return(NhlFATAL);
}
}
}
/*
* Get argument # 4
*/
qvp = (void*)NclGetArgValue(
4,
8,
&ndims_qvp,
dsizes_qvp,
NULL,
NULL,
&type_qvp,
DONT_CARE);
/*
* Check dimension sizes.
*/
if(ndims_qvp != ndims_pres) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: The qvp and pres arrays must have the same dimensionality");
return(NhlFATAL);
}
else {
for(i = 0; i < ndims_pres; i++) {
if(dsizes_qvp[i] != dsizes_pres[i]) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: The qvp and pres arrays must have the same dimensionality");
return(NhlFATAL);
}
}
}
/*
* Get argument # 5
*/
ght = (void*)NclGetArgValue(
5,
8,
&ndims_ght,
dsizes_ght,
NULL,
NULL,
&type_ght,
DONT_CARE);
/*
* Check dimension sizes.
*/
if(ndims_ght != ndims_pres) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: The ght and pres arrays must have the same dimensionality");
return(NhlFATAL);
}
else {
for(i = 0; i < ndims_pres; i++) {
if(dsizes_ght[i] != dsizes_pres[i]) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: The ght and pres arrays must have the same dimensionality");
return(NhlFATAL);
}
}
}
/*
* Get argument # 6
*/
ter = (void*)NclGetArgValue(
6,
8,
&ndims_ter,
dsizes_ter,
NULL,
NULL,
&type_ter,
DONT_CARE);
/*
* Check dimension sizes for ter. It can either be 2D, or one fewer
* dimensions than pres.
*/
if(ndims_ter != 2 && ndims_ter != (ndims_pres-1)) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: ter must either be a 2D array dimensioned south_north x west_east or it must have the same dimensionality as the pres array, minus the level dimension");
return(NhlFATAL);
}
if(ndims_ter == 2) {
if(dsizes_ter[0] != nlat || dsizes_ter[1] != nlon) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: The dimensions of ter must be south_north x west_east");
return(NhlFATAL);
}
}
else {
for(i = 0; i < ndims_pres-3; i++) {
if(dsizes_ter[i] != dsizes_pres[i]) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: ter must either be a 2D array dimensioned south_north x west_east or it must have the same dimensionality as the pres array, minus the level dimension");
return(NhlFATAL);
}
}
}
/*
* Get argument # 7
*/
haveqci = (int*)NclGetArgValue(
7,
8,
NULL,
NULL,
NULL,
NULL,
NULL,
DONT_CARE);
/*
* Calculate size of leftmost dimensions.
*/
size_leftmost = 1;
for(i = 0; i < ndims_pres-3; i++) size_leftmost *= dsizes_pres[i];
/*
* Allocate space for coercing input arrays. If any of the input
* is already double, then we don't need to allocate space for
* temporary arrays, because we'll just change the pointer into
* the void array appropriately.
*/
/*
* Allocate space for tmp_pres.
*/
nlatlon = nlat * nlon;
nlevlatlon = nlev * nlatlon;
if(type_pres != NCL_double) {
tmp_pres = (double *)calloc(nlevlatlon,sizeof(double));
if(tmp_pres == NULL) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: Unable to allocate memory for coercing pressure array to double");
return(NhlFATAL);
}
}
/*
* Allocate space for tmp_tk.
*/
if(type_tk != NCL_double) {
tmp_tk = (double *)calloc(nlevlatlon,sizeof(double));
if(tmp_tk == NULL) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: Unable to allocate memory for coercing tk array to double");
return(NhlFATAL);
}
}
/*
* Allocate space for tmp_qci.
*/
if(type_qci != NCL_double) {
tmp_qci = (double *)calloc(nlevlatlon,sizeof(double));
if(tmp_qci == NULL) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: Unable to allocate memory for coercing qci array to double");
return(NhlFATAL);
}
}
/*
* Allocate space for tmp_qcw.
*/
if(type_qcw != NCL_double) {
tmp_qcw = (double *)calloc(nlevlatlon,sizeof(double));
if(tmp_qcw == NULL) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: Unable to allocate memory for coercing qcw array to double");
return(NhlFATAL);
}
}
/*
* Allocate space for tmp_qvp.
*/
if(type_qvp != NCL_double) {
tmp_qvp = (double *)calloc(nlevlatlon,sizeof(double));
if(tmp_qvp == NULL) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: Unable to allocate memory for coercing qvp array to double");
return(NhlFATAL);
}
}
/*
* Allocate space for tmp_ght.
*/
if(type_ght != NCL_double) {
tmp_ght = (double *)calloc(nlevlatlon,sizeof(double));
if(tmp_ght == NULL) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: Unable to allocate memory for coercing ght array to double");
return(NhlFATAL);
}
}
/*
* Coerce ter to double, if necessary.
*/
if(ndims_ter == 2) {
tmp_ter = coerce_input_double(ter,type_ter,nlatlon,0,NULL,NULL);
if(tmp_ter == NULL) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: Unable to allocate memory for coercing ter array to double");
return(NhlFATAL);
}
}
else {
/*
* Allocate space for tmp_ter.
*/
if(type_ter != NCL_double) {
tmp_ter = (double *)calloc(nlatlon,sizeof(double));
if(tmp_ter == NULL) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: Unable to allocate memory for coercing ter array to double");
return(NhlFATAL);
}
}
}
/*
* The output type defaults to float, unless one or more input
* arrays are double.
*/
if(type_pres == NCL_double || type_tk == NCL_double ||
type_qci == NCL_double || type_qcw == NCL_double ||
type_qvp == NCL_double || type_ght == NCL_double ||
type_ter == NCL_double) {
type_ctt = NCL_double;
type_obj_ctt = nclTypedoubleClass;
}
else {
type_ctt = NCL_float;
type_obj_ctt = nclTypefloatClass;
}
/*
* Allocate space for output array.
*/
size_output = size_leftmost * nlatlon;
if(type_ctt != NCL_double) {
ctt = (void *)calloc(size_output, sizeof(float));
tmp_ctt = (double *)calloc(nlatlon,sizeof(double));
if(ctt == NULL || tmp_ctt == NULL) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: Unable to allocate memory for temporary output array");
return(NhlFATAL);
}
}
else {
ctt = (void *)calloc(size_output, sizeof(double));
if(ctt == NULL) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: Unable to allocate memory for output array");
return(NhlFATAL);
}
}
/*
* Allocate space for output dimension sizes and set them.
*/
ndims_ctt = ndims_pres-1;
dsizes_ctt = (ng_size_t*)calloc(ndims_ctt,sizeof(ng_size_t));
if( dsizes_ctt == NULL ) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: Unable to allocate memory for holding dimension sizes");
return(NhlFATAL);
}
for(i = 0; i < ndims_ctt-2; i++) dsizes_ctt[i] = dsizes_pres[i];
dsizes_ctt[ndims_ctt-2] = nlat;
dsizes_ctt[ndims_ctt-1] = nlon;
/*
* Get dimension info to see if we have named dimensions.
* Using "ght" here, because it is more likely than "pres"
* to have metadata attached to it.
*
* This will be used for return variable.
*/
dim_info_ght = get_wrf_dim_info(5,8,ndims_ght,dsizes_ght);
if(dim_info_ght != NULL) {
dim_info = malloc(sizeof(NclDimRec)*ndims_ctt);
if(dim_info == NULL) {
NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_ctt: Unable to allocate memory for holding dimension information");
return(NhlFATAL);
}
for(i = 0; i < ndims_ght-3; i++) {
dim_info[i] = dim_info_ght[i];
}
dim_info[ndims_ctt-1] = dim_info_ght[ndims_ght-1];
dim_info[ndims_ctt-2] = dim_info_ght[ndims_ght-2];
}
/*
* Loop across leftmost dimensions and call the Fortran routine for each
* subsection of the input arrays.
*/
index_pres = index_ter = index_ctt = 0;
for(i = 0; i < size_leftmost; i++) {
/*
* Coerce subsection of pres (tmp_pres) to double if necessary.
*/
if(type_pres != NCL_double) {
coerce_subset_input_double(pres,tmp_pres,index_pres,
type_pres,nlevlatlon,0,NULL,NULL);
}
else {
tmp_pres = &((double*)pres)[index_pres];
}
/*
* Coerce subsection of tk (tmp_tk) to double if necessary.
*/
if(type_tk != NCL_double) {
coerce_subset_input_double(tk,tmp_tk,index_pres,type_tk,
nlevlatlon,0,NULL,NULL);
}
else {
tmp_tk = &((double*)tk)[index_pres];
}
/*
* Coerce subsection of qci (tmp_qci) to double if necessary.
*/
if(type_qci != NCL_double) {
coerce_subset_input_double(qci,tmp_qci,index_pres,type_qci,
nlevlatlon,0,NULL,NULL);
}
else {
tmp_qci = &((double*)qci)[index_pres];
}
/*
* Coerce subsection of qcw (tmp_qcw) to double if necessary.
*/
if(type_qcw != NCL_double) {
coerce_subset_input_double(qcw,tmp_qcw,index_pres,type_qcw,
nlevlatlon,0,NULL,NULL);
}
else {
tmp_qcw = &((double*)qcw)[index_pres];
}
/*
* Coerce subsection of qvp (tmp_qvp) to double if necessary.
*/
if(type_qvp != NCL_double) {
coerce_subset_input_double(qvp,tmp_qvp,index_pres,type_qvp,
nlevlatlon,0,NULL,NULL);
}
else {
tmp_qvp = &((double*)qvp)[index_pres];
}
/*
* Coerce subsection of ght (tmp_ght) to double if necessary.
*/
if(type_ght != NCL_double) {
coerce_subset_input_double(ght,tmp_ght,index_pres,type_ght,
nlevlatlon,0,NULL,NULL);
}
else {
tmp_ght = &((double*)ght)[index_pres];
}
/*
* Coerce subsection of ter (tmp_ter) to double if necessary.
*/
if(ndims_ter != 2) {
if(type_ter != NCL_double) {
coerce_subset_input_double(ter,tmp_ter,index_ter,type_ter,
nlatlon,0,NULL,NULL);
}
else {
tmp_ter = &((double*)ter)[index_ter];
}
}
/*
* Point temporary output array to void output array if appropriate.
*/
if(type_ctt == NCL_double) {
tmp_ctt = &((double*)ctt)[index_ctt];
}
/*
* Call the Fortran routine.
*/
NGCALLF(wrfcttcalc,WRFCTTCALC)(tmp_pres, tmp_tk, tmp_qci, tmp_qcw,
tmp_qvp, tmp_ght, tmp_ter, tmp_ctt,
haveqci,&inlev, &inlat, &inlon);
/*
* Coerce output back to float if necessary.
*/
if(type_ctt == NCL_float) {
coerce_output_float_only(ctt,tmp_ctt,nlatlon,
index_ctt);
}
index_pres += nlevlatlon;
index_ctt += nlatlon;
if(ndims_ter != 2) {
index_ter += nlatlon;
}
}
/*
* Free unneeded memory.
*/
if(type_pres != NCL_double) NclFree(tmp_pres);
if(type_tk != NCL_double) NclFree(tmp_tk);
if(type_qci != NCL_double) NclFree(tmp_qci);
if(type_qcw != NCL_double) NclFree(tmp_qcw);
if(type_qvp != NCL_double) NclFree(tmp_qvp);
if(type_ght != NCL_double) NclFree(tmp_ght);
if(type_ter != NCL_double) NclFree(tmp_ter);
if(type_ctt != NCL_double) NclFree(tmp_ctt);
/*
* Set up some attributes ("description" and "units") to return.
*/
cdescription = (char *)calloc(22,sizeof(char));
cunits = (char *)calloc(2,sizeof(char));
strcpy(cdescription,"Cloud Top Temperature");
strcpy(cunits,"K");
description = (NclQuark*)NclMalloc(sizeof(NclQuark));
units = (NclQuark*)NclMalloc(sizeof(NclQuark));
*description = NrmStringToQuark(cdescription);
*units = NrmStringToQuark(cunits);
free(cdescription);
free(cunits);
/*
* Set up return value.
*/
return_md = _NclCreateVal(
NULL,
NULL,
Ncl_MultiDValData,
0,
(void*)ctt,
NULL,
ndims_ctt,
dsizes_ctt,
TEMPORARY,
NULL,
type_obj_ctt
);
/*
* Set up attributes to return.
*/
att_id = _NclAttCreate(NULL,NULL,Ncl_Att,0,NULL);
dsizes[0] = 1;
att_md = _NclCreateVal(
NULL,
NULL,
Ncl_MultiDValData,
0,
(void*)description,
NULL,
1,
dsizes,
TEMPORARY,
NULL,
(NclObjClass)nclTypestringClass
);
_NclAddAtt(
att_id,
"description",
att_md,
NULL
);
att_md = _NclCreateVal(
NULL,
NULL,
Ncl_MultiDValData,
0,
(void*)units,
NULL,
1,
dsizes,
TEMPORARY,
NULL,
(NclObjClass)nclTypestringClass
);
_NclAddAtt(
att_id,
"units",
att_md,
NULL
);
tmp_var = _NclVarCreate(
NULL,
NULL,
Ncl_Var,
0,
NULL,
return_md,
dim_info,
att_id,
NULL,
RETURNVAR,
NULL,
TEMPORARY
);
if(dim_info != NULL) NclFree(dim_info);
NclFree(dim_info_ght);
/*
* Return output grid and attributes to NCL.
*/
return_data.kind = NclStk_VAR;
return_data.u.data_var = tmp_var;
_NclPlaceReturn(return_data);
return(NhlNOERROR);
}

117
ncl_reference/wrf_fctt.f
Unescape View File

@ -1,117 +0,0 @@
C NCLFORTSTART
subroutine wrfcttcalc(prs,tk,qci,qcw,qvp,ght,ter,ctt,
& haveqci,nz,ns,ew)
implicit none
integer nz,ns,ew,haveqci
double precision ght(ew,ns,nz)
double precision prs(ew,ns,nz),tk(ew,ns,nz)
double precision qci(ew,ns,nz),qcw(ew,ns,nz)
double precision qvp(ew,ns,nz)
double precision ctt(ew,ns),ter(ew,ns)
c double precision znfac(nz)
C NCLEND
c
c
integer i,j,k,mjx,miy,mkzh,ripk,wrfout
double precision vt,rgas,grav,opdepthu,opdepthd,dp
double precision ratmix,eps,arg1,arg2,agl_hgt,ussalr
double precision abscoefi,abscoef,fac,prsctt,celkel
c double precision ght(ew,ns,nz),stuff(ew,ns)
double precision pf(ns,ew,nz),p1,p2
c
c
mjx = ew
miy = ns
mkzh = nz
eps = 0.622d0
ussalr = .0065d0 ! deg C per m
rgas = 287.04d0 !J/K/kg
grav = 9.81d0
abscoefi = .272d0 ! cloud ice absorption coefficient in m^2/g
abscoef =.145d0 ! cloud water absorption coefficient in m^2/g
celkel = 273.15d0
wrfout = 1
cCalculate the surface pressure
do j=1,ew
do i=1,ns
ratmix = .001d0*qvp(j,i,1)
arg1 = eps + ratmix
arg2 = eps*(1.+ratmix)
vt = tk(j,i,1) * arg1/arg2 !Virtual temperature
agl_hgt = ght(j,i,nz) - ter(j,i)
arg1 = -grav/(rgas*ussalr)
pf(i,j,nz) = prs(j,i,1)*
& (vt/(vt+ussalr*(agl_hgt)))**(arg1)
enddo
enddo
c
do j=1,ew
do i=1,ns
do k=1,nz-1
ripk = nz-k+1
pf(i,j,k)=.5d0*(prs(j,i,ripk)+prs(j,i,ripk-1))
enddo
enddo
enddo
do 190 j=1,ew
do 190 i=1,ns
opdepthd=0.d0
k=0
c
c Integrate downward from model top, calculating path at full
c model vertical levels.
c
20 opdepthu=opdepthd
k=k+1
ripk = nz-k+1
if (k.eq.1) then
dp=200.d0*(pf(i,j,1)-prs(j,i,nz)) ! should be in Pa
else
dp=100.d0*(pf(i,j,k)-pf(i,j,k-1)) ! should be in Pa
endif
if (haveqci .eq. 0) then
if (tk(i,j,k).lt.celkel) then
c Note: abscoefi is m**2/g, qcw is g/kg,
c so no convrsion needed
opdepthd=opdepthu+abscoefi*qcw(j,i,k)*dp/grav
else
opdepthd=opdepthu+abscoef*qcw(j,i,k)*dp/grav
endif
else
opdepthd=opdepthd+(abscoef*qcw(j,i,ripk)+
& abscoefi*qci(j,i,ripk))*dp/grav
endif
if (opdepthd.lt.1..and.k.lt.nz) then
goto 20
elseif (opdepthd.lt.1..and.k.eq.nz) then
prsctt=prs(j,i,1)
else
fac=(1.-opdepthu)/(opdepthd-opdepthu)
prsctt=pf(i,j,k-1)+fac*(pf(i,j,k)-pf(i,j,k-1))
prsctt=min(prs(j,i,1),max(prs(j,i,nz),prsctt))
endif
do 30 k=2,nz
ripk = nz-k+1
p1 = prs(j,i,ripk+1)
p2 = prs(j,i,ripk)
if (prsctt .ge. p1 .and. prsctt .le .p2) then
fac=(prsctt-p1)/(p2-p1)
arg1 = fac*(tk(j,i,ripk)-tk(j,i,ripk+1))-celkel
ctt(j,i) = tk(j,i,ripk+1)+ arg1
goto 40
endif
30 continue
40 continue
190 continue
return
end

1917
ncl_reference/wrf_fddaobs_in.F
View File

File diff suppressed because it is too large Load Diff

380
ncl_reference/wrf_map_fix.ncl
Unescape View File

@ -1,380 +0,0 @@
undef("set_mp_wrf_map_resources")
function set_mp_wrf_map_resources(in_file[1]:file,opt_args[1]:logical)
begin
;
opts = opt_args ; Make a copy of the resource list
; Set some resources depending on what kind of map projection is
; chosen.
;
; MAP_PROJ = 0 : "CylindricalEquidistant"
; MAP_PROJ = 1 : "LambertConformal"
; MAP_PROJ = 2 : "Stereographic"
; MAP_PROJ = 3 : "Mercator"
; MAP_PROJ = 6 : "Lat/Lon"
if(isatt(in_file,"MAP_PROJ"))
; CylindricalEquidistant
if(in_file@MAP_PROJ .eq. 0)
projection = "CylindricalEquidistant"
opts@mpProjection = projection
opts@mpGridSpacingF = 45
opts@mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF", 0.0)
if(isatt(in_file,"STAND_LON"))
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",in_file@STAND_LON)
else
if(isatt(in_file,"CEN_LON"))
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",in_file@CEN_LON)
else
print("ERROR: Found neither STAND_LON or CEN_LON in file")
end if
end if
end if
; LambertConformal projection
if(in_file@MAP_PROJ .eq. 1)
projection = "LambertConformal"
opts@mpProjection = projection
opts@mpLambertParallel1F = get_res_value_keep(opts, "mpLambertParallel1F",in_file@TRUELAT1)
opts@mpLambertParallel2F = get_res_value_keep(opts, "mpLambertParallel2F",in_file@TRUELAT2)
if(isatt(in_file,"STAND_LON"))
opts@mpLambertMeridianF = get_res_value_keep(opts, "mpLambertMeridianF",in_file@STAND_LON)
else
if(isatt(in_file,"CEN_LON"))
opts@mpLambertMeridianF = get_res_value_keep(opts, "mpLambertMeridianF",in_file@CEN_LON)
else
print("ERROR: Found neither STAND_LON or CEN_LON in file")
end if
end if
end if
; Stereographic projection
if(in_file@MAP_PROJ .eq. 2)
projection = "Stereographic"
opts@mpProjection = projection
opts@mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF", in_file@CEN_LAT)
if(isatt(in_file,"STAND_LON"))
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",in_file@STAND_LON)
else
if(isatt(in_file,"CEN_LON"))
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",in_file@CEN_LON)
else
print("ERROR: Found neither STAND_LON or CEN_LON in file")
end if
end if
end if
; Mercator projection
if(in_file@MAP_PROJ .eq. 3)
projection = "Mercator"
opts@mpProjection = projection
opts@mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF", 0.0)
if(isatt(in_file,"STAND_LON"))
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",in_file@STAND_LON)
else
if(isatt(in_file,"CEN_LON"))
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",in_file@CEN_LON)
else
print("ERROR: Found neither STAND_LON or CEN_LON in file")
end if
end if
end if
; global WRF CylindricalEquidistant
if(in_file@MAP_PROJ .eq. 6)
projection = "CylindricalEquidistant"
opts@mpProjection = projection
opts@mpGridSpacingF = 45
if (isatt(in_file,"POLE_LON") .and. isatt(in_file,"POLE_LAT") .and. isatt(in_file,"STAND_LON")) then
if (in_file@POLE_LON .eq. 0 .and. in_file@POLE_LAT .eq. 90) then
; not rotated
opts@mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF", 0.0)
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",180 - in_file@STAND_LON)
else
; rotated
southern = False ; default to northern hemisphere
if (in_file@POLE_LON .eq. 0.0) then
southern = True
else if (in_file@POLE_LON .ne. 180) then
if (isatt(in_file,"CEN_LAT") .and. in_file@CEN_LAT .lt. 0.0) then
southern = True ; probably but not necessarily true -- no way to tell for sure
end if
end if
end if
if (.not. southern) then
opts@mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF", 90.0 - in_file@POLE_LAT)
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF", -in_file@STAND_LON)
else
opts@mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF", in_file@POLE_LAT - 90)
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF", 180 - in_file@STAND_LON)
end if
end if
else if (isatt(in_file,"ref_lon") .and. isatt(in_file,"ref_lat")) then
;; this is definitely true for NMM grids but unlikely for ARW grids especially if ref_x and ref_y are set
opts@mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF", in_file@REF_LAT)
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF", in_file@REF_LON)
else if (isatt(in_file,"cen_lat") .and. isatt(in_file,"cen_lon")) then
;; these usually specifiy the center of the coarse domain --- not necessarily the center of the projection
opts@mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF",in_file@CEN_LAT)
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",in_file@CEN_LON)
else
;; default values for global grid
opts@mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF", 0.0)
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF", 180.0)
end if
end if
end if
end if
end if
return(opts) ; Return.
end
undef("wrf_map_resources")
function wrf_map_resources(in_file[1]:file,map_args[1]:logical)
local lat, lon, x1, x2, y1, y2, dims, ii, jj, southern
begin
;
; This function sets resources for a WRF map plot, basing the projection on
; the MAP_PROJ attribute in the given file. It's intended to be callable
; by users who need to set mpXXXX resources for other plotting scripts.
;
; Set some resources depending on what kind of map projection is
; chosen.
;
; MAP_PROJ = 0 : "CylindricalEquidistant"
; MAP_PROJ = 1 : "LambertConformal"
; MAP_PROJ = 2 : "Stereographic"
; MAP_PROJ = 3 : "Mercator"
; MAP_PROJ = 6 : "Lat/Lon"
if(isatt(in_file,"MAP_PROJ"))
; CylindricalEquidistant
if(in_file@MAP_PROJ .eq. 0)
map_args@mpProjection = "CylindricalEquidistant"
map_args@mpGridSpacingF = 45
map_args@mpCenterLatF = get_res_value_keep(map_args, "mpCenterLatF", 0.0)
if(isatt(in_file,"STAND_LON"))
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF",in_file@STAND_LON)
else
if(isatt(in_file,"CEN_LON"))
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF",in_file@CEN_LON)
else
print("ERROR: Found neither STAND_LON or CEN_LON in file")
end if
end if
end if
; LambertConformal projection
if(in_file@MAP_PROJ .eq. 1)
map_args@mpProjection = "LambertConformal"
map_args@mpLambertParallel1F = get_res_value_keep(map_args, "mpLambertParallel1F",in_file@TRUELAT1)
map_args@mpLambertParallel2F = get_res_value_keep(map_args, "mpLambertParallel2F",in_file@TRUELAT2)
if(isatt(in_file,"STAND_LON"))
map_args@mpLambertMeridianF = get_res_value_keep(map_args, "mpLambertMeridianF",in_file@STAND_LON)
else
if(isatt(in_file,"CEN_LON"))
map_args@mpLambertMeridianF = get_res_value_keep(map_args, "mpLambertMeridianF",in_file@CEN_LON)
else
print("ERROR: Found neither STAND_LON or CEN_LON in file")
end if
end if
end if
; Stereographic projection
if(in_file@MAP_PROJ .eq. 2)
map_args@mpProjection = "Stereographic"
map_args@mpCenterLatF = get_res_value_keep(map_args, "mpCenterLatF", in_file@CEN_LAT)
if(isatt(in_file,"STAND_LON"))
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF",in_file@STAND_LON)
else
if(isatt(in_file,"CEN_LON"))
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF",in_file@CEN_LON)
else
print("ERROR: Found neither STAND_LON or CEN_LON in file")
end if
end if
end if
; Mercator projection
if(in_file@MAP_PROJ .eq. 3)
map_args@mpProjection = "Mercator"
map_args@mpCenterLatF = get_res_value_keep(map_args, "mpCenterLatF", 0.0)
if(isatt(in_file,"STAND_LON"))
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF",in_file@STAND_LON)
else
if(isatt(in_file,"CEN_LON"))
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF",in_file@CEN_LON)
else
print("ERROR: Found neither STAND_LON or CEN_LON in file")
end if
end if
end if
; global WRF CylindricalEquidistant
if(in_file@MAP_PROJ .eq. 6)
print ("YES, THIS WORKED")
projection = "CylindricalEquidistant"
map_args@mpProjection = projection
map_args@mpGridSpacingF = 45
;; according to the docs if POLE_LON is 0 then the projection center is in the southern hemisphere
;; if POLE_LON is 180 the projection center is in the northern hemisphere
;; otherwise you can't tell for sure -- CEN_LAT does not have to be the projection center but hopefully
;; it is in the same hemisphere. The same is true for REF_LAT except that if REF_Y is specified REF_LAT might
;; be in a corner or somewhere else and therefore it is even less reliable
;;
if (isatt(in_file,"POLE_LON") .and. isatt(in_file,"POLE_LAT") .and. isatt(in_file,"STAND_LON")) then
if (in_file@POLE_LON .eq. 0 .and. in_file@POLE_LAT .eq. 90) then
; not rotated
map_args@mpCenterLatF = get_res_value_keep(map_args, "mpCenterLatF", 0.0)
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF",180 - in_file@STAND_LON)
else
; rotated
southern = False ; default to northern hemisphere
if (in_file@POLE_LON .eq. 0.0) then
southern = True
else if (in_file@POLE_LON .ne. 180) then
if (isatt(in_file,"CEN_LAT") .and. in_file@CEN_LAT .lt. 0.0) then
southern = True ; probably but not necessarily true -- no way to tell for sure
end if
end if
end if
if (.not. southern) then
map_args@mpCenterLatF = get_res_value_keep(map_args, "mpCenterLatF", 90.0 - in_file@POLE_LAT)
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF", -in_file@STAND_LON)
else
map_args@mpCenterLatF = get_res_value_keep(map_args, "mpCenterLatF", in_file@POLE_LAT - 90)
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF", 180 - in_file@STAND_LON)
end if
end if
else if (isatt(in_file,"ref_lon") .and. isatt(in_file,"ref_lat")) then
;; this is definitely true for NMM grids but unlikely for ARW grids especially if ref_x and ref_y are set
map_args@mpCenterLatF = get_res_value_keep(map_args, "mpCenterLatF", in_file@REF_LAT)
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF", in_file@REF_LON)
else if (isatt(in_file,"cen_lat") .and. isatt(in_file,"cen_lon")) then
;; these usually specifiy the center of the coarse domain --- not necessarily the center of the projection
map_args@mpCenterLatF = get_res_value_keep(map_args, "mpCenterLatF",in_file@CEN_LAT)
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF",in_file@CEN_LON)
else
;; default values for global grid
map_args@mpCenterLatF = get_res_value_keep(map_args, "mpCenterLatF", 0.0)
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF", 180.0)
end if
end if
end if
end if
else
return(map_args)
end if
map_args@mpNestTime = get_res_value_keep(map_args, "mpNestTime",0)
if(isfilevar(in_file,"XLAT"))
lat = in_file->XLAT(map_args@mpNestTime,:,:)
lon = in_file->XLONG(map_args@mpNestTime,:,:)
else
lat = in_file->XLAT_M(map_args@mpNestTime,:,:)
lon = in_file->XLONG_M(map_args@mpNestTime,:,:)
end if
dims = dimsizes(lat)
do ii = 0, dims(0)-1
do jj = 0, dims(1)-1
if ( lon(ii,jj) .lt. 0.0) then
lon(ii,jj) = lon(ii,jj) + 360.
end if
end do
end do
map_args@start_lat = lat(0,0)
map_args@start_lon = lon(0,0)
map_args@end_lat = lat(dims(0)-1,dims(1)-1)
map_args@end_lon = lon(dims(0)-1,dims(1)-1)
; end_lon must be greater than start_lon, or errors are thrown
if (map_args@end_lon .le. map_args@start_lon) then
map_args@end_lon = map_args@end_lon + 360.0
end if
; Set some resources common to all map projections.
map_args = set_mp_resources(map_args)
if ( isatt(map_args,"ZoomIn") .and. map_args@ZoomIn ) then
y1 = 0
x1 = 0
y2 = dims(0)-1
x2 = dims(1)-1
if ( isatt(map_args,"Ystart") ) then
y1 = map_args@Ystart
delete(map_args@Ystart)
end if
if ( isatt(map_args,"Xstart") ) then
x1 = map_args@Xstart
delete(map_args@Xstart)
end if
if ( isatt(map_args,"Yend") ) then
if ( map_args@Yend .le. y2 ) then
y2 = map_args@Yend
end if
delete(map_args@Yend)
end if
if ( isatt(map_args,"Xend") ) then
if ( map_args@Xend .le. x2 ) then
x2 = map_args@Xend
end if
delete(map_args@Xend)
end if
map_args@mpLeftCornerLatF = lat(y1,x1)
map_args@mpLeftCornerLonF = lon(y1,x1)
map_args@mpRightCornerLatF = lat(y2,x2)
map_args@mpRightCornerLonF = lon(y2,x2)
if ( map_args@mpRightCornerLonF .lt. 0.0 ) then
map_args@mpRightCornerLonF = map_args@mpRightCornerLonF + 360.0
end if
if ( map_args@mpRightCornerLonF .le. map_args@mpRightCornerLonF ) then
map_args@mpRightCornerLonF = map_args@mpRightCornerLonF + 360.0
end if
delete(map_args@ZoomIn)
end if
return(map_args)
end

109
ncl_reference/wrf_pvo.f
Unescape View File

@ -1,109 +0,0 @@
c--------------------------------------------------------
C NCLFORTSTART
SUBROUTINE DCOMPUTEPV(PV,U,V,THETA,PRS,MSFU,MSFV,MSFT,COR,DX,DY,
+ NX,NY,NZ,NXP1,NYP1)
IMPLICIT NONE
INTEGER NX,NY,NZ,NXP1,NYP1
DOUBLE PRECISION U(NXP1,NY,NZ),V(NX,NYP1,NZ),PRS(NX,NY,NZ)
DOUBLE PRECISION THETA(NX,NY,NZ),PV(NX,NY,NZ)
DOUBLE PRECISION MSFU(NXP1,NY),MSFV(NX,NYP1),MSFT(NX,NY)
DOUBLE PRECISION COR(NX,NY)
DOUBLE PRECISION DX,DY
C NCLEND
INTEGER KP1,KM1,JP1,JM1,IP1,IM1,I,J,K
DOUBLE PRECISION DSY,DSX,DP,DUDY,DVDX,DUDP,DVDP,DTHDP,AVORT
DOUBLE PRECISION DTHDX,DTHDY,MM
c print*,'nx,ny,nz,nxp1,nyp1'
c print*,nx,ny,nz,nxp1,nyp1
DO K = 1,NZ
KP1 = MIN(K+1,NZ)
KM1 = MAX(K-1,1)
DO J = 1,NY
JP1 = MIN(J+1,NY)
JM1 = MAX(J-1,1)
DO I = 1,NX
IP1 = MIN(I+1,NX)
IM1 = MAX(I-1,1)
c print *,jp1,jm1,ip1,im1
DSX = (IP1-IM1)*DX
DSY = (JP1-JM1)*DY
MM = MSFT(I,J)*MSFT(I,J)
c print *,j,i,u(i,jp1,k),msfu(i,jp1),u(i,jp1,k)/msfu(i,jp1)
DUDY = 0.5D0* (U(I,JP1,K)/MSFU(I,JP1)+
+ U(I+1,JP1,K)/MSFU(I+1,JP1)-
+ U(I,JM1,K)/MSFU(I,JM1)-
+ U(I+1,JM1,K)/MSFU(I+1,JM1))/DSY*MM
DVDX = 0.5D0* (V(IP1,J,K)/MSFV(IP1,J)+
+ V(IP1,J+1,K)/MSFV(IP1,J+1)-
+ V(IM1,J,K)/MSFV(IM1,J)-
+ V(IM1,J+1,K)/MSFV(IM1,J+1))/DSX*MM
AVORT = DVDX - DUDY + COR(I,J)
DP = PRS(I,J,KP1) - PRS(I,J,KM1)
DUDP = 0.5D0* (U(I,J,KP1)+U(I+1,J,KP1)-U(I,J,KM1)-
+ U(I+1,J,KM1))/DP
DVDP = 0.5D0* (V(I,J,KP1)+V(I,J+1,KP1)-V(I,J,KM1)-
+ V(I,J+1,KM1))/DP
DTHDP = (THETA(I,J,KP1)-THETA(I,J,KM1))/DP
DTHDX = (THETA(IP1,J,K)-THETA(IM1,J,K))/DSX*MSFT(I,J)
DTHDY = (THETA(I,JP1,K)-THETA(I,JM1,K))/DSY*MSFT(I,J)
PV(I,J,K) = -9.81D0* (DTHDP*AVORT-DVDP*DTHDX+
+ DUDP*DTHDY)*10000.D0
c if(i.eq.300 .and. j.eq.300) then
c print*,'avort,dudp,dvdp,dthdp,dthdx,dthdy,pv'
c print*,avort,dudp,dvdp,dthdp,dthdx,dthdy,pv(i,j,k)
c endif
PV(I,J,K) = PV(I,J,K)*1.D2
END DO
END DO
END DO
RETURN
END
c--------------------------------------------------------
C NCLFORTSTART
SUBROUTINE DCOMPUTEABSVORT(AV,U,V,MSFU,MSFV,MSFT,COR,DX,DY,NX,NY,
+ NZ,NXP1,NYP1)
IMPLICIT NONE
INTEGER NX,NY,NZ,NXP1,NYP1
DOUBLE PRECISION U(NXP1,NY,NZ),V(NX,NYP1,NZ)
DOUBLE PRECISION AV(NX,NY,NZ)
DOUBLE PRECISION MSFU(NXP1,NY),MSFV(NX,NYP1),MSFT(NX,NY)
DOUBLE PRECISION COR(NX,NY)
DOUBLE PRECISION DX,DY
C NCLEND
INTEGER KP1,KM1,JP1,JM1,IP1,IM1,I,J,K
DOUBLE PRECISION DSY,DSX,DP,DUDY,DVDX,DUDP,DVDP,DTHDP,AVORT
DOUBLE PRECISION DTHDX,DTHDY,MM
c print*,'nx,ny,nz,nxp1,nyp1'
c print*,nx,ny,nz,nxp1,nyp1
DO K = 1,NZ
DO J = 1,NY
JP1 = MIN(J+1,NY)
JM1 = MAX(J-1,1)
DO I = 1,NX
IP1 = MIN(I+1,NX)
IM1 = MAX(I-1,1)
c print *,jp1,jm1,ip1,im1
DSX = (IP1-IM1)*DX
DSY = (JP1-JM1)*DY
MM = MSFT(I,J)*MSFT(I,J)
c print *,j,i,u(i,jp1,k),msfu(i,jp1),u(i,jp1,k)/msfu(i,jp1)
DUDY = 0.5D0* (U(I,JP1,K)/MSFU(I,JP1)+
+ U(I+1,JP1,K)/MSFU(I+1,JP1)-
+ U(I,JM1,K)/MSFU(I,JM1)-
+ U(I+1,JM1,K)/MSFU(I+1,JM1))/DSY*MM
DVDX = 0.5D0* (V(IP1,J,K)/MSFV(IP1,J)+
+ V(IP1,J+1,K)/MSFV(IP1,J+1)-
+ V(IM1,J,K)/MSFV(IM1,J)-
+ V(IM1,J+1,K)/MSFV(IM1,J+1))/DSX*MM
AVORT = DVDX - DUDY + COR(I,J)
AV(I,J,K) = AVORT*1.D5
END DO
END DO
END DO
RETURN
END

100
ncl_reference/wrf_relhl.f
Unescape View File

@ -1,100 +0,0 @@
C ***************************************************************
C * Storm Relative Helicity (SRH) is a measure of the *
C * streamwise vorticity within the inflow environment of a *
C * convective storm. It is calculated by multiplying the *
C * storm-relative inflow velocity vector (Vh-C) by the *
C * streamwise vorticity (Zh) and integrating this quantity *
C * over the inflow depth (lowest 1-3 km layers above ground *
C * level). It describes the extent to which corkscrew-like *
C * motion occurs (similar to the spiraling motion of an *
C * American football). SRH corresponds to the transfer of *
C * vorticity from the environment to an air parcel in *
C * convective motion and is used to predict the potential *
C * for tornadic development (cyclonic updraft rotation) in *
C * right-moving supercells. *
C * *
C * There is no clear threshold value for SRH when forecasting *
C * supercells, since the formation of supercells appears to be *
C * related more strongly to the deeper layer vertical shear. *
C * Larger values of 0-3-km SRH (greater than 250 m**2/s**2) *
C * and 0-1-km SRH (greater than 100 m**2/s**2), suggest an *
C * increased threat of tornadoes with supercells. For SRH, *
C * larger values are generally better, but there are no clear *
C * "boundaries" between non-tornadic and significant tornadic *
C * supercells. *
C * *
C * SRH < 100 (lowest 1 km): cutoff value *
C * SRH = 150-299: supercells possible with weak tornadoes *
C * SRH = 300-499: very favorable to supercell development and *
C * strong tornadoes *
C * SRH > 450 : violent tornadoes *
C ***************************************************************
C NCLFORTSTART
subroutine dcalrelhl(u, v, ght, ter, top, sreh, miy, mjx, mkzh)
implicit none
integer miy, mjx, mkzh
double precision u(miy,mjx,mkzh), v(miy,mjx,mkzh),
& ght(miy,mjx,mkzh),top,ter(miy,mjx),
& sreh(miy,mjx)
C NCLEND
C
C This helicity code was provided by Dr. Craig Mattocks, and
C verified by Cindy Bruyere to produce results equivalent to
C those generated by RIP4. (The code came from RIP4?)
C
double precision pi, dtr, dpr
double precision dh, sdh, su, sv, ua, va, asp, adr, bsp, bdr
double precision cu, cv, x, sum
integer i, j, k, k10, k3, ktop
parameter (pi=3.14159265d0, dtr=pi/180.d0, dpr=180.d0/pi)
do 15 j = 1, mjx-1
do 15 i = 1, miy-1
sdh = 0.d0
su = 0.d0
sv = 0.d0
k3 = 0
k10 = 0
ktop = 0
do 6 k = mkzh, 2, -1
if (((ght(i,j,k) - ter(i,j)) .gt. 10000.d0) .and.
& (k10 .eq. 0)) then
k10 = k
go to 8
endif
if (((ght(i,j,k) - ter(i,j)) .gt. top) .and.
& (ktop .eq. 0)) ktop = k
if (((ght(i,j,k) - ter(i,j)) .gt. 3000.d0) .and.
& (k3 .eq. 0)) k3 = k
6 continue
8 continue
if (k10 .eq. 0) k10=2
do k = k3, k10, -1
dh = ght(i,j,k-1) - ght(i,j,k)
sdh = sdh + dh
su = su + 0.5d0*dh*(u(i,j,k-1)+u(i,j,k))
sv = sv + 0.5d0*dh*(v(i,j,k-1)+v(i,j,k))
enddo
ua = su / sdh
va = sv / sdh
asp = sqrt(ua*ua + va*va)
if (ua .eq. 0.d0 .and. va .eq. 0.d0) then
adr = 0.d0
else
adr = dpr * (pi + atan2(ua,va))
endif
bsp = 0.75d0 * asp
bdr = adr + 30.d0
if (bdr .gt. 360.d0) bdr = bdr-360.d0
cu = -bsp * sin(bdr*dtr)
cv = -bsp * cos(bdr*dtr)
sum = 0.d0
do 12 k = mkzh-1, ktop, -1
x = ((u(i,j,k)-cu) * (v(i,j,k)-v(i,j,k+1))) -
& ((v(i,j,k)-cv) * (u(i,j,k)-u(i,j,k+1)))
sum = sum + x
12 continue
sreh(i,j) = -sum
15 continue
end

264
ncl_reference/wrf_rip_phys_routines.f
Unescape View File

@ -1,264 +0,0 @@
c======================================================================
c
c !IROUTINE: WETBULBCALC -- Calculate wet bulb temperature (C)
c
c !DESCRIPTION:
c
c Calculates wet bulb temperature in C, given pressure in
c temperature in K and mixing ratio in kg/kg.
c
c !INPUT:
c nx - index for x dimension
c ny - index for y dimension
c nz - index for z dimension
c prs - pressure (mb)
c tmk - temperature (K)
c qvp - water vapor mixing ratio (kg/kg)
c
c !OUTPUT:
c twb - Wet bulb temperature (C)
c
c !ASSUMPTIONS:
c
c !REVISION HISTORY:
c 2009-March - Mark T. Stoelinga - from RIP4.5
c 2010-August - J. Schramm
c 2014-March - A. Jaye - modified to run with NCL and ARW wrf output
c
c !INTERFACE:
c ------------------------------------------------------------------
C NCLFORTSTART
subroutine wetbulbcalc(prs,tmk,qvp,twb,nx,ny,nz,psafile)
implicit none
integer nx, ny, nz
double precision prs(nz,ny,nx)
double precision tmk(nz,ny,nx)
double precision qvp(nz,ny,nx)
double precision twb(nz,ny,nx)
character*(*) psafile
C NCLEND
integer i,j,k
integer jtch,jt,ipch,ip
double precision q, t, p, e, tlcl, eth
double precision fracip,fracip2,fracjt,fracjt2
double precision PSADITHTE(150),PSADIPRS(150),PSADITMK(150,150)
double precision tonpsadiabat
double precision eps,tlclc1,tlclc2,tlclc3,tlclc4,gamma
double precision gammamd,thtecon1,thtecon2,thtecon3,celkel
double precision rgas,rgasmd,cp,cpmd
c
c Before looping, set lookup table for getting temperature on
c a pseudoadiabat.
c
CALL DLOOKUP_TABLE(PSADITHTE,PSADIPRS,PSADITMK,psafile)
c Define constants
rgas=287.04 !J/K/kg
rgasmd=.608 ! rgas_moist=rgas*(1.+rgasmd*qvp)
cp=1004. ! J/K/kg Note: not using Bolton's value of 1005.7
cpmd=.887 ! cp_moist=cp*(1.+cpmd*qvp)
eps=0.622
tlclc1=2840.
tlclc2=3.5
tlclc3=4.805
tlclc4=55.
gamma=rgas/cp
gammamd=rgasmd-cpmd ! gamma_moist=gamma*(1.+gammamd*qvp)
thtecon1=3376. ! K
thtecon2=2.54
thtecon3=.81
celkel=273.15
DO k=1,nx
DO j=1,ny
DO i=1,nz
q=dmax1(qvp(i,j,k),1.d-15)
t=tmk(i,j,k)
p=prs(i,j,k)/100.
e=q*p/(eps+q)
tlcl=tlclc1/(dlog(t**tlclc2/e)-tlclc3)+tlclc4
eth=t*(1000./p)**(gamma*(1.+gammamd*q))*
& exp((thtecon1/tlcl-thtecon2)*q*(1.+thtecon3*q))
c
c Now we need to find the temperature (in K) on a moist adiabat
c (specified by eth in K) given pressure in hPa. It uses a
c lookup table, with data that was generated by the Bolton (1980)
c formula for theta_e.
c
c First check if pressure is less than min pressure in lookup table.
c If it is, assume parcel is so dry that the given theta-e value can
c be interpretted as theta, and get temperature from the simple dry
c theta formula.
c
if (p.le.psadiprs(150)) then
tonpsadiabat=eth*(p/1000.)**gamma
else
c
c Otherwise, look for the given thte/prs point in the lookup table.
c
do jtch=1,150-1
if (eth.ge.psadithte(jtch).and.eth.lt.
& psadithte(jtch+1)) then
jt=jtch
goto 213
endif
enddo
jt=-1
213 continue
do ipch=1,150-1
if (p.le.psadiprs(ipch).and.p.gt.psadiprs(ipch+1)) then
ip=ipch
goto 215
endif
enddo
ip=-1
215 continue
if (jt.eq.-1.or.ip.eq.-1) then
print*,
& 'Outside of lookup table bounds. prs,thte=',p,eth
stop
endif
fracjt=(eth-psadithte(jt))/(psadithte(jt+1)-psadithte(jt))
fracjt2=1.-fracjt
fracip=(psadiprs(ip)-p)/(psadiprs(ip)-psadiprs(ip+1))
fracip2=1.-fracip
if (psaditmk(ip,jt).gt.1e9.or.psaditmk(ip+1,jt).gt.1e9.or.
& psaditmk(ip,jt+1).gt.1e9.or.
& psaditmk(ip+1,jt+1).gt.1e9) then
print*,
& 'Tried to access missing tmperature in lookup table.'
print*,
& 'Prs and Thte probably unreasonable. prs,thte='
& ,p,eth
stop
endif
tonpsadiabat=fracip2*fracjt2*psaditmk(ip ,jt )+
& fracip *fracjt2*psaditmk(ip+1,jt )+
& fracip2*fracjt *psaditmk(ip ,jt+1)+
& fracip *fracjt *psaditmk(ip+1,jt+1)
endif
twb(i,j,k)=tonpsadiabat
ENDDO
ENDDO
ENDDO
c
return
end
c======================================================================
c
c !IROUTINE: omgcalc -- Calculate omega (dp/dt)
c
c !DESCRIPTION:
c
c Calculate approximate omega, based on vertical velocity w (dz/dt).
c It is approximate because it cannot take into account the vertical
c motion of pressure surfaces.
c
c !INPUT:
c mx - index for x dimension
c my - index for y dimension
c mx - index for vertical dimension
c qvp - water vapor mixing ratio (kg/kg)
c tmk - temperature (K)
c www - vertical velocity (m/s)
c prs - pressure (Pa)
c
c !OUTPUT:
c omg - omega (Pa/sec)
c
c !ASSUMPTIONS:
c
c !REVISION HISTORY:
c 2009-March - Mark T. Stoelinga - from RIP4.5
c 2010-August - J. Schramm
c 2014-March - A. Jaye - modified to run with NCL and ARW wrf output
c
c ------------------------------------------------------------------
c NCLFORTSTART
subroutine omgcalc(qvp,tmk,www,prs,omg,mx,my,mz)
implicit none
integer mx, my, mz
double precision qvp(mz,my,mx)
double precision tmk(mz,my,mx)
double precision www(mz,my,mx)
double precision prs(mz,my,mx)
double precision omg(mz,my,mx)
c NCLEND
c Local variables
integer i, j, k
double precision grav,rgas,eps
c
c Constants
c
grav=9.81 ! m/s**2
rgas=287.04 !J/K/kg
eps=0.622
do k=1,mx
do j=1,my
do i=1,mz
omg(i,j,k)=-grav*prs(i,j,k)/
& (rgas*((tmk(i,j,k)*(eps+qvp(i,j,k)))/
& (eps*(1.+qvp(i,j,k)))))*www(i,j,k)
enddo
enddo
enddo
c
return
end
c======================================================================
c
c !IROUTINE: VIRTUAL_TEMP -- Calculate virtual temperature (K)
c
c !DESCRIPTION:
c
c Calculates virtual temperature in K, given temperature
c in K and mixing ratio in kg/kg.
c
c !INPUT:
c NX - index for x dimension
c NY - index for y dimension
c NZ - index for z dimension
c RATMIX - water vapor mixing ratio (kg/kg)
c TEMP - temperature (K)
c
c !OUTPUT:
c TV - Virtual temperature (K)
c
c !ASSUMPTIONS:
c
c !REVISION HISTORY:
c 2009-March - Mark T. Stoelinga - from RIP4.5
c 2010-August - J. Schramm
c 2014-March - A. Jaye - modified to run with NCL and ARW wrf output
c
c ------------------------------------------------------------------
C NCLFORTSTART
SUBROUTINE VIRTUAL_TEMP(TEMP,RATMIX,TV,NX,NY,NZ)
IMPLICIT NONE
INTEGER NX,NY,NZ
DOUBLE PRECISION TEMP(NZ,NY,NX)
DOUBLE PRECISION RATMIX(NZ,NY,NX)
DOUBLE PRECISION TV(NZ,NY,NX)
C NCLEND
INTEGER I,J,K
DOUBLE PRECISION EPS
EPS = 0.622D0
DO K=1,NX
DO J=1,NY
DO I=1,NZ
TV(I,J,K) = TEMP(I,J,K)* (EPS+RATMIX(I,J,K))/
& (EPS* (1.D0+RATMIX(I,J,K)))
ENDDO
ENDDO
ENDDO
RETURN
END

771
ncl_reference/wrf_user.f
Unescape View File

@ -1,771 +0,0 @@
C NCLFORTSTART
SUBROUTINE DCOMPUTEPI(PI,PRESSURE,NX,NY,NZ)
IMPLICIT NONE
INTEGER NX,NY,NZ
DOUBLE PRECISION PI(NX,NY,NZ)
DOUBLE PRECISION PRESSURE(NX,NY,NZ)
C NCLEND
INTEGER I,J,K
DOUBLE PRECISION P1000MB,R_D,CP
PARAMETER (P1000MB=100000.D0,R_D=287.D0,CP=7.D0*R_D/2.D0)
DO K = 1,NZ
DO J = 1,NY
DO I = 1,NX
PI(I,J,K) = (PRESSURE(I,J,K)/P1000MB)** (R_D/CP)
END DO
END DO
END DO
END
C NCLFORTSTART
SUBROUTINE DCOMPUTETK(TK,PRESSURE,THETA,NX)
IMPLICIT NONE
INTEGER NX
DOUBLE PRECISION PI
DOUBLE PRECISION PRESSURE(NX)
DOUBLE PRECISION THETA(NX)
DOUBLE PRECISION TK(NX)
C NCLEND
INTEGER I
DOUBLE PRECISION P1000MB,R_D,CP
PARAMETER (P1000MB=100000.D0,R_D=287.D0,CP=7.D0*R_D/2.D0)
DO I = 1,NX
PI = (PRESSURE(I)/P1000MB)** (R_D/CP)
TK(I) = PI*THETA(I)
END DO
END
C NCLFORTSTART
SUBROUTINE DINTERP3DZ(V3D,V2D,Z,LOC,NX,NY,NZ,VMSG)
IMPLICIT NONE
INTEGER NX,NY,NZ
DOUBLE PRECISION V3D(NX,NY,NZ),V2D(NX,NY)
DOUBLE PRECISION Z(NX,NY,NZ)
DOUBLE PRECISION LOC
DOUBLE PRECISION VMSG
C NCLEND
INTEGER I,J,KP,IP,IM
LOGICAL INTERP
DOUBLE PRECISION HEIGHT,W1,W2
HEIGHT = LOC
c does vertical coordinate increase or decrease with increasing k?
c set offset appropriately
IP = 0
IM = 1
IF (Z(1,1,1).GT.Z(1,1,NZ)) THEN
IP = 1
IM = 0
END IF
DO I = 1,NX
DO J = 1,NY
C Initialize to missing. Was initially hard-coded to -999999.
V2D(I,J) = VMSG
INTERP = .false.
KP = NZ
DO WHILE ((.NOT.INTERP) .AND. (KP.GE.2))
IF (((Z(I,J,KP-IM).LE.HEIGHT).AND. (Z(I,J,
+ KP-IP).GT.HEIGHT))) THEN
W2 = (HEIGHT-Z(I,J,KP-IM))/
+ (Z(I,J,KP-IP)-Z(I,J,KP-IM))
W1 = 1.D0 - W2
V2D(I,J) = W1*V3D(I,J,KP-IM) + W2*V3D(I,J,KP-IP)
INTERP = .true.
END IF
KP = KP - 1
END DO
END DO
END DO
RETURN
END
C NCLFORTSTART
SUBROUTINE DZSTAG(ZNEW,NX,NY,NZ,Z,NXZ,NYZ,NZZ,TERRAIN)
IMPLICIT NONE
INTEGER NX,NY,NZ,NXZ,NYZ,NZZ
DOUBLE PRECISION ZNEW(NX,NY,NZ),Z(NXZ,NYZ,NZZ)
DOUBLE PRECISION TERRAIN(NXZ,NYZ)
C NCLEND
INTEGER I,J,K,II,IM1,JJ,JM1
c check for u, v, or w (x,y,or z) staggering
c
c for x and y stag, avg z to x, y, point
c
IF (NX.GT.NXZ) THEN
DO K = 1,NZ
DO J = 1,NY
DO I = 1,NX
II = MIN0(I,NXZ)
IM1 = MAX0(I-1,1)
ZNEW(I,J,K) = 0.5D0* (Z(II,J,K)+Z(IM1,J,K))
END DO
END DO
END DO
ELSE IF (NY.GT.NYZ) THEN
DO K = 1,NZ
DO J = 1,NY
JJ = MIN0(J,NYZ)
JM1 = MAX0(J-1,1)
DO I = 1,NX
ZNEW(I,J,K) = 0.5D0* (Z(I,JJ,K)+Z(I,JM1,K))
END DO
END DO
END DO
c
c w (z) staggering
c
ELSE IF (NZ.GT.NZZ) THEN
DO J = 1,NY
DO I = 1,NX
ZNEW(I,J,1) = TERRAIN(I,J)
END DO
END DO
DO K = 2,NZ
DO J = 1,NY
DO I = 1,NX
ZNEW(I,J,K) = ZNEW(I,J,K-1) +
+ 2.D0* (Z(I,J,K-1)-ZNEW(I,J,K-1))
END DO
END DO
END DO
END IF
RETURN
END
C NCLFORTSTART
SUBROUTINE DINTERP2DXY(V3D,V2D,XY,NX,NY,NZ,NXY)
IMPLICIT NONE
INTEGER NX,NY,NZ,NXY
DOUBLE PRECISION V3D(NX,NY,NZ),V2D(NXY,NZ)
DOUBLE PRECISION XY(2,NXY)
C NCLEND
INTEGER I,J,K,IJ
DOUBLE PRECISION W11,W12,W21,W22,WX,WY
DO IJ = 1,NXY
I = MAX0(1,MIN0(NX-1,INT(XY(1,IJ)+1)))
J = MAX0(1,MIN0(NY-1,INT(XY(2,IJ)+1)))
WX = DBLE(I+1) - (XY(1,IJ)+1)
WY = DBLE(J+1) - (XY(2,IJ)+1)
W11 = WX*WY
W21 = (1.D0-WX)*WY
W12 = WX* (1.D0-WY)
W22 = (1.D0-WX)* (1.D0-WY)
DO K = 1,NZ
V2D(IJ,K) = W11*V3D(I,J,K) + W21*V3D(I+1,J,K) +
+ W12*V3D(I,J+1,K) + W22*V3D(I+1,J+1,K)
END DO
END DO
RETURN
END
C NCLFORTSTART
SUBROUTINE DINTERP1D(V_IN,V_OUT,Z_IN,Z_OUT,NZ_IN,NZ_OUT,VMSG)
IMPLICIT NONE
INTEGER NZ_IN,NZ_OUT
DOUBLE PRECISION V_IN(NZ_IN),Z_IN(NZ_IN)
DOUBLE PRECISION V_OUT(NZ_OUT),Z_OUT(NZ_OUT)
DOUBLE PRECISION VMSG
C NCLEND
INTEGER KP,K,IM,IP
LOGICAL INTERP
DOUBLE PRECISION HEIGHT,W1,W2
c does vertical coordinate increase of decrease with increasing k?
c set offset appropriately
IP = 0
IM = 1
IF (Z_IN(1).GT.Z_IN(NZ_IN)) THEN
IP = 1
IM = 0
END IF
DO K = 1,NZ_OUT
V_OUT(K) = VMSG
INTERP = .false.
KP = NZ_IN
HEIGHT = Z_OUT(K)
DO WHILE ((.NOT.INTERP) .AND. (KP.GE.2))
IF (((Z_IN(KP-IM).LE.HEIGHT).AND.
+ (Z_IN(KP-IP).GT.HEIGHT))) THEN
W2 = (HEIGHT-Z_IN(KP-IM))/ (Z_IN(KP-IP)-Z_IN(KP-IM))
W1 = 1.D0 - W2
V_OUT(K) = W1*V_IN(KP-IM) + W2*V_IN(KP-IP)
INTERP = .true.
END IF
KP = KP - 1
END DO
END DO
RETURN
END
c---------------------------------------------
c Bill,
c This routine assumes
c index order is (i,j,k)
c wrf staggering
C
c units: pressure (Pa), temperature(K), height (m), mixing ratio
c (kg kg{-1}) availability of 3d p, t, and qv; 2d terrain; 1d
c half-level zeta string
c output units of SLP are Pa, but you should divide that by 100 for the
c weather weenies.
c virtual effects are included
c
C NCLFORTSTART
SUBROUTINE DCOMPUTESEAPRS(NX,NY,NZ,Z,T,P,Q,SEA_LEVEL_PRESSURE,
+ T_SEA_LEVEL,T_SURF,LEVEL)
IMPLICIT NONE
c Estimate sea level pressure.
INTEGER NX,NY,NZ
DOUBLE PRECISION Z(NX,NY,NZ)
DOUBLE PRECISION T(NX,NY,NZ),P(NX,NY,NZ),Q(NX,NY,NZ)
c The output is the 2d sea level pressure.
DOUBLE PRECISION SEA_LEVEL_PRESSURE(NX,NY)
INTEGER LEVEL(NX,NY)
DOUBLE PRECISION T_SURF(NX,NY),T_SEA_LEVEL(NX,NY)
C NCLEND
c Some required physical constants:
DOUBLE PRECISION R,G,GAMMA
PARAMETER (R=287.04D0,G=9.81D0,GAMMA=0.0065D0)
c Specific constants for assumptions made in this routine:
DOUBLE PRECISION TC,PCONST
PARAMETER (TC=273.16D0+17.5D0,PCONST=10000)
LOGICAL RIDICULOUS_MM5_TEST
PARAMETER (RIDICULOUS_MM5_TEST=.TRUE.)
c PARAMETER (ridiculous_mm5_test = .false.)
c Local variables:
INTEGER I,J,K
INTEGER KLO,KHI
DOUBLE PRECISION PLO,PHI,TLO,THI,ZLO,ZHI
DOUBLE PRECISION P_AT_PCONST,T_AT_PCONST,Z_AT_PCONST
DOUBLE PRECISION Z_HALF_LOWEST
LOGICAL L1,L2,L3,FOUND
C
c Find least zeta level that is PCONST Pa above the surface. We
c later use this level to extrapolate a surface pressure and
c temperature, which is supposed to reduce the effect of the diurnal
c heating cycle in the pressure field.
DO J = 1,NY
DO I = 1,NX
LEVEL(I,J) = -1
K = 1
FOUND = .false.
DO WHILE ((.NOT.FOUND) .AND. (K.LE.NZ))
IF (P(I,J,K).LT.P(I,J,1)-PCONST) THEN
LEVEL(I,J) = K
FOUND = .true.
END IF
K = K + 1
END DO
IF (LEVEL(I,J).EQ.-1) THEN
PRINT '(A,I4,A)','Troubles finding level ',
+ NINT(PCONST)/100,' above ground.'
PRINT '(A,I4,A,I4,A)','Problems first occur at (',I,
+ ',',J,')'
PRINT '(A,F6.1,A)','Surface pressure = ',P(I,J,1)/100,
+ ' hPa.'
STOP 'Error_in_finding_100_hPa_up'
END IF
END DO
END DO
c Get temperature PCONST Pa above surface. Use this to extrapolate
c the temperature at the surface and down to sea level.
DO J = 1,NY
DO I = 1,NX
KLO = MAX(LEVEL(I,J)-1,1)
KHI = MIN(KLO+1,NZ-1)
IF (KLO.EQ.KHI) THEN
PRINT '(A)','Trapping levels are weird.'
PRINT '(A,I3,A,I3,A)','klo = ',KLO,', khi = ',KHI,
+ ': and they should not be equal.'
STOP 'Error_trapping_levels'
END IF
PLO = P(I,J,KLO)
PHI = P(I,J,KHI)
TLO = T(I,J,KLO)* (1.D0+0.608D0*Q(I,J,KLO))
THI = T(I,J,KHI)* (1.D0+0.608D0*Q(I,J,KHI))
c zlo = zetahalf(klo)/ztop*(ztop-terrain(i,j))+terrain(i,j)
c zhi = zetahalf(khi)/ztop*(ztop-terrain(i,j))+terrain(i,j)
ZLO = Z(I,J,KLO)
ZHI = Z(I,J,KHI)
P_AT_PCONST = P(I,J,1) - PCONST
T_AT_PCONST = THI - (THI-TLO)*LOG(P_AT_PCONST/PHI)*
+ LOG(PLO/PHI)
Z_AT_PCONST = ZHI - (ZHI-ZLO)*LOG(P_AT_PCONST/PHI)*
+ LOG(PLO/PHI)
T_SURF(I,J) = T_AT_PCONST* (P(I,J,1)/P_AT_PCONST)**
+ (GAMMA*R/G)
T_SEA_LEVEL(I,J) = T_AT_PCONST + GAMMA*Z_AT_PCONST
END DO
END DO
C
c If we follow a traditional computation, there is a correction to the
c sea level temperature if both the surface and sea level
c temperatures are *too* hot.
IF (RIDICULOUS_MM5_TEST) THEN
DO J = 1,NY
DO I = 1,NX
L1 = T_SEA_LEVEL(I,J) .LT. TC
L2 = T_SURF(I,J) .LE. TC
L3 = .NOT. L1
IF (L2 .AND. L3) THEN
T_SEA_LEVEL(I,J) = TC
ELSE
T_SEA_LEVEL(I,J) = TC -
+ 0.005D0* (T_SURF(I,J)-TC)**2
END IF
END DO
END DO
END IF
c The grand finale: ta da!
DO J = 1,NY
DO I = 1,NX
c z_half_lowest=zetahalf(1)/ztop*(ztop-terrain(i,j))+terrain(i,j)
Z_HALF_LOWEST = Z(I,J,1)
C Convert to hPa in this step, by multiplying by 0.01. The original
C Fortran routine didn't do this, but the NCL script that called it
C did, so we moved it here.
SEA_LEVEL_PRESSURE(I,J) = 0.01 * (P(I,J,1)*
+ EXP((2.D0*G*Z_HALF_LOWEST)/
+ (R* (T_SEA_LEVEL(I,J)+T_SURF(I,
+ J)))))
END DO
END DO
c print *,'sea pres input at weird location i=20,j=1,k=1'
c print *,'t=',t(20,1,1),t(20,2,1),t(20,3,1)
c print *,'z=',z(20,1,1),z(20,2,1),z(20,3,1)
c print *,'p=',p(20,1,1),p(20,2,1),p(20,3,1)
c print *,'slp=',sea_level_pressure(20,1),
c * sea_level_pressure(20,2),sea_level_pressure(20,3)
END
c---------------------------------------------------
C
C Double precision version. If you make a change here, you
C must make the same change below to filter2d.
C
C NCLFORTSTART
SUBROUTINE DFILTER2D(A,B,NX,NY,IT,MISSING)
IMPLICIT NONE
c Estimate sea level pressure.
INTEGER NX,NY,IT
DOUBLE PRECISION A(NX,NY),B(NX,NY),MISSING
C NCLEND
DOUBLE PRECISION COEF
PARAMETER (COEF=0.25D0)
INTEGER I,J,ITER
DO ITER = 1,IT
DO J = 1,NY
DO I = 1,NX
B(I,J) = A(I,J)
END DO
END DO
DO J = 2,NY - 1
DO I = 1,NX
IF ( B(I,J-1).EQ.MISSING .OR. B(I,J).EQ.MISSING .OR.
+ B(I,J+1).EQ.MISSING ) THEN
A(I,J) = A(I,J)
ELSE
A(I,J) = A(I,J) + COEF* (B(I,J-1)-2*B(I,J)+B(I,J+1))
END IF
END DO
END DO
DO J = 1,NY
DO I = 2,NX - 1
IF ( B(I-1,J).EQ.MISSING .OR. B(I,J).EQ.MISSING .OR.
+ B(I+1,J).EQ.MISSING ) THEN
A(I,J) = A(I,J)
ELSE
A(I,J) = A(I,J) + COEF* (B(I-1,J)-2*B(I,J)+B(I+1,J))
END IF
END DO
END DO
c do j=1,ny
c do i=1,nx
c b(i,j) = a(i,j)
c enddo
c enddo
c do j=2,ny-1
c do i=1,nx
c a(i,j) = a(i,j) - .99*coef*(b(i,j-1)-2*b(i,j)+b(i,j+1))
c enddo
c enddo
c do j=1,ny
c do i=2,nx-1
c a(i,j) = a(i,j) - .99*coef*(b(i-1,j)-2*b(i,j)+b(i+1,j))
c enddo
c enddo
END DO
RETURN
END
C
C Single precision version. If you make a change here, you
C must make the same change above to dfilter2d.
C
C NCLFORTSTART
SUBROUTINE filter2d( a, b, nx , ny , it, missing)
IMPLICIT NONE
c Estimate sea level pressure.
INTEGER nx , ny, it
REAL a(nx,ny),b(nx,ny), missing
C NCLEND
REAL coef
parameter( coef = 0.25)
INTEGER i,j,iter
do iter=1, it
do j=1,ny
do i=1,nx
b(i,j) = a(i,j)
enddo
enddo
do j=2,ny-1
do i=1,nx
if ( b(i,j-1).eq.missing .or. b(i,j).eq.missing .or.
+ b(i,j+1).eq.missing ) then
a(i,j) = a(i,j)
else
a(i,j) = a(i,j) + coef*(b(i,j-1)-2*b(i,j)+b(i,j+1))
end if
enddo
enddo
do j=1,ny
do i=2,nx-1
if ( b(i-1,j).eq.missing .or. b(i,j).eq.missing .or.
+ b(i+1,j).eq.missing ) then
a(i,j) = a(i,j)
else
a(i,j) = a(i,j) + coef*(b(i-1,j)-2*b(i,j)+b(i+1,j))
end if
enddo
enddo
c do j=1,ny
c do i=1,nx
c b(i,j) = a(i,j)
c enddo
c enddo
c do j=2,ny-1
c do i=1,nx
c a(i,j) = a(i,j) - .99*coef*(b(i,j-1)-2*b(i,j)+b(i,j+1))
c enddo
c enddo
c do j=1,ny
c do i=2,nx-1
c a(i,j) = a(i,j) - .99*coef*(b(i-1,j)-2*b(i,j)+b(i+1,j))
c enddo
c enddo
enddo
return
end
c---------------------------------------------------------
C NCLFORTSTART
SUBROUTINE DCOMPUTERH(QV,P,T,RH,NX)
IMPLICIT NONE
INTEGER NX
DOUBLE PRECISION QV(NX),P(NX),T(NX),RH(NX)
C NCLEND
DOUBLE PRECISION SVP1,SVP2,SVP3,SVPT0
PARAMETER (SVP1=0.6112D0,SVP2=17.67D0,SVP3=29.65D0,SVPT0=273.15D0)
INTEGER I
DOUBLE PRECISION QVS,ES,PRESSURE,TEMPERATURE
DOUBLE PRECISION EP_2,R_D,R_V
PARAMETER (R_D=287.D0,R_V=461.6D0,EP_2=R_D/R_V)
DOUBLE PRECISION EP_3
PARAMETER (EP_3=0.622D0)
DO I = 1,NX
PRESSURE = P(I)
TEMPERATURE = T(I)
c es = 1000.*svp1*
ES = 10.D0*SVP1*EXP(SVP2* (TEMPERATURE-SVPT0)/
+ (TEMPERATURE-SVP3))
c qvs = ep_2*es/(pressure-es)
QVS = EP_3*ES/ (0.01D0*PRESSURE- (1.D0-EP_3)*ES)
c rh = 100*amax1(1., qv(i)/qvs)
c rh(i) = 100.*qv(i)/qvs
RH(I) = 100.D0*DMAX1(DMIN1(QV(I)/QVS,1.0D0),0.0D0)
END DO
RETURN
END
c----------------------------------------------
C NCLFORTSTART
SUBROUTINE DGETIJLATLONG(LAT_ARRAY,LONG_ARRAY,LAT,LONGITUDE,
+ II,JJ,NX,NY,IMSG)
IMPLICIT NONE
INTEGER NX,NY,II,JJ,IMSG
DOUBLE PRECISION LAT_ARRAY(NX,NY),LONG_ARRAY(NX,NY)
DOUBLE PRECISION LAT,LONGITUDE
C NCLEND
DOUBLE PRECISION LONGD,LATD
INTEGER I,J
DOUBLE PRECISION IR,JR
DOUBLE PRECISION DIST_MIN,DIST
C Init to missing. Was hard-coded to -999 initially.
IR = IMSG
JR = IMSG
DIST_MIN = 1.D+20
DO J = 1,NY
DO I = 1,NX
LATD = (LAT_ARRAY(I,J)-LAT)**2
LONGD = (LONG_ARRAY(I,J)-LONGITUDE)**2
C LONGD = DMIN1((LONG_ARRAY(I,J)-LONGITUDE)**2,
C + (LONG_ARRAY(I,J)+LONGITUDE)**2)
DIST = SQRT(LATD+LONGD)
IF (DIST_MIN.GT.DIST) THEN
DIST_MIN = DIST
IR = DBLE(I)
JR = DBLE(J)
END IF
END DO
END DO
C
C The original version of this routine returned IR and JR. But, then
C the NCL script that called this routine was converting IR and JR
C to integer, so why not just return II and JJ?
C
C Also, I'm subtracing 1 here, because it will be returned to NCL
C script which has 0-based indexing.
C
IF(IR.ne.IMSG.and.JR.ne.IMSG) then
II = NINT(IR)-1
JJ = NINT(JR)-1
ELSE
II = IMSG
JJ = IMSG
END IF
c we will just return the nearest point at present
RETURN
END
C NCLFORTSTART
SUBROUTINE DCOMPUTEUVMET(U,V,UVMET,LONGCA,LONGCB,FLONG,FLAT,
+ CEN_LONG,CONE,RPD,NX,NY,NXP1,NYP1,
+ ISTAG,IS_MSG_VAL,UMSG,VMSG,UVMETMSG)
IMPLICIT NONE
C ISTAG should be 0 if the U,V grids are not staggered.
C That is, NY = NYP1 and NX = NXP1.
INTEGER NX,NY,NXP1,NYP1,ISTAG
LOGICAL IS_MSG_VAL
DOUBLE PRECISION U(NXP1,NY),V(NX,NYP1)
DOUBLE PRECISION UVMET(NX,NY,2)
DOUBLE PRECISION FLONG(NX,NY),FLAT(NX,NY)
DOUBLE PRECISION LONGCB(NX,NY),LONGCA(NX,NY)
DOUBLE PRECISION CEN_LONG,CONE,RPD
DOUBLE PRECISION UMSG,VMSG,UVMETMSG
C NCLEND
INTEGER I,J
DOUBLE PRECISION UK,VK
c WRITE (6,FMT=*) ' in compute_uvmet ',NX,NY,NXP1,NYP1,ISTAG
DO J = 1,NY
DO I = 1,NX
LONGCA(I,J) = FLONG(I,J) - CEN_LONG
IF (LONGCA(I,J).GT.180.D0) THEN
LONGCA(I,J) = LONGCA(I,J) - 360.D0
END IF
IF (LONGCA(I,J).LT.-180.D0) THEN
LONGCA(I,J) = LONGCA(I,J) + 360.D0
END IF
IF (FLAT(I,J).LT.0.D0) THEN
LONGCB(I,J) = -LONGCA(I,J)*CONE*RPD
ELSE
LONGCB(I,J) = LONGCA(I,J)*CONE*RPD
END IF
LONGCA(I,J) = COS(LONGCB(I,J))
LONGCB(I,J) = SIN(LONGCB(I,J))
END DO
END DO
c WRITE (6,FMT=*) ' computing velocities '
DO J = 1,NY
DO I = 1,NX
IF (ISTAG.EQ.1) THEN
IF (IS_MSG_VAL.AND.(U(I,J).EQ.UMSG.OR.
+ V(I,J).EQ.VMSG.OR.
+ U(I+1,J).EQ.UMSG.OR.
+ V(I,J+1).EQ.VMSG)) THEN
UVMET(I,J,1) = UVMETMSG
UVMET(I,J,2) = UVMETMSG
ELSE
UK = 0.5D0* (U(I,J)+U(I+1,J))
VK = 0.5D0* (V(I,J)+V(I,J+1))
UVMET(I,J,1) = VK*LONGCB(I,J) + UK*LONGCA(I,J)
UVMET(I,J,2) = VK*LONGCA(I,J) - UK*LONGCB(I,J)
END IF
ELSE
IF (IS_MSG_VAL.AND.(U(I,J).EQ.UMSG.OR.
+ V(I,J).EQ.VMSG)) THEN
UVMET(I,J,1) = UVMETMSG
UVMET(I,J,2) = UVMETMSG
ELSE
UK = U(I,J)
VK = V(I,J)
UVMET(I,J,1) = VK*LONGCB(I,J) + UK*LONGCA(I,J)
UVMET(I,J,2) = VK*LONGCA(I,J) - UK*LONGCB(I,J)
END IF
END IF
END DO
END DO
RETURN
END
C NCLFORTSTART
C
C This was originally a routine that took 2D input arrays. Since
C the NCL C wrapper routine can handle multiple dimensions, it's
C not necessary to have anything bigger than 1D here.
C
SUBROUTINE DCOMPUTETD(TD,PRESSURE,QV_IN,NX)
IMPLICIT NONE
INTEGER NX
DOUBLE PRECISION PRESSURE(NX)
DOUBLE PRECISION QV_IN(NX)
DOUBLE PRECISION TD(NX)
C NCLEND
DOUBLE PRECISION QV,TDC
INTEGER I
DO I = 1,NX
QV = DMAX1(QV_IN(I),0.D0)
c vapor pressure
TDC = QV*PRESSURE(I)/ (.622D0+QV)
c avoid problems near zero
TDC = DMAX1(TDC,0.001D0)
TD(I) = (243.5D0*LOG(TDC)-440.8D0)/ (19.48D0-LOG(TDC))
END DO
RETURN
END
C NCLFORTSTART
SUBROUTINE DCOMPUTEICLW(ICLW,PRESSURE,QC_IN,NX,NY,NZ)
IMPLICIT NONE
INTEGER NX,NY,NZ
DOUBLE PRECISION PRESSURE(NX,NY,NZ)
DOUBLE PRECISION QC_IN(NX,NY,NZ)
DOUBLE PRECISION ICLW(NX,NY)
DOUBLE PRECISION QCLW,DP,GG
C NCLEND
INTEGER I,J,K
GG = 1000.D0/9.8D0
DO J = 1,NY
DO I = 1,NX
ICLW(I,J) = 0.D0
END DO
END DO
DO J = 3,NY - 2
DO I = 3,NX - 2
DO K = 1,NZ
QCLW = DMAX1(QC_IN(I,J,K),0.D0)
IF (K.EQ.1) THEN
DP = (PRESSURE(I,J,K-1)-PRESSURE(I,J,K))
ELSE IF (K.EQ.NZ) THEN
DP = (PRESSURE(I,J,K)-PRESSURE(I,J,K+1))
ELSE
DP = (PRESSURE(I,J,K-1)-PRESSURE(I,J,K+1))/2.D0
END IF
ICLW(I,J) = ICLW(I,J) + QCLW*DP*GG
END DO
END DO
END DO
RETURN
END

209
ncl_reference/wrf_user_dbz.f
Unescape View File

@ -1,209 +0,0 @@
C NCLFORTSTART
SUBROUTINE CALCDBZ(DBZ,PRS,TMK,QVP,QRA,QSN,QGR,WEDIM,SNDIM,BTDIM,
+ SN0,IVARINT,ILIQSKIN)
c
c This routine computes equivalent reflectivity factor (in dBZ) at
c each model grid point. In calculating Ze, the RIP algorithm makes
c assumptions consistent with those made in an early version
c (ca. 1996) of the bulk mixed-phase microphysical scheme in the MM5
c model (i.e., the scheme known as "Resiner-2"). For each species:
c
c 1. Particles are assumed to be spheres of constant density. The
c densities of rain drops, snow particles, and graupel particles are
c taken to be rho_r = rho_l = 1000 kg m^-3, rho_s = 100 kg m^-3, and
c rho_g = 400 kg m^-3, respectively. (l refers to the density of
c liquid water.)
c
c 2. The size distribution (in terms of the actual diameter of the
c particles, rather than the melted diameter or the equivalent solid
c ice sphere diameter) is assumed to follow an exponential
c distribution of the form N(D) = N_0 * exp( lambda*D ).
c
c 3. If ivarint=0, the intercept parameters are assumed constant
c (as in early Reisner-2), with values of 8x10^6, 2x10^7,
c and 4x10^6 m^-4, for rain, snow, and graupel, respectively.
c If ivarint=1, variable intercept parameters are used, as
c calculated in Thompson, Rasmussen, and Manning (2004, Monthly
c Weather Review, Vol. 132, No. 2, pp. 519-542.)
c
c 4. If iliqskin=1, frozen particles that are at a temperature above
c freezing are assumed to scatter as a liquid particle.
c
c More information on the derivation of simulated reflectivity in
c RIP can be found in Stoelinga (2005, unpublished write-up).
c Contact Mark Stoelinga (stoeling@atmos.washington.edu) for a copy.
c
c Arguments
INTEGER WEDIM,SNDIM,BTDIM
INTEGER SN0,IVARINT,ILIQSKIN
DOUBLE PRECISION DBZ(WEDIM,SNDIM,BTDIM)
DOUBLE PRECISION PRS(WEDIM,SNDIM,BTDIM)
DOUBLE PRECISION TMK(WEDIM,SNDIM,BTDIM)
DOUBLE PRECISION QVP(WEDIM,SNDIM,BTDIM)
DOUBLE PRECISION QRA(WEDIM,SNDIM,BTDIM)
DOUBLE PRECISION QSN(WEDIM,SNDIM,BTDIM)
DOUBLE PRECISION QGR(WEDIM,SNDIM,BTDIM)
C NCLEND
c Local Variables
INTEGER I,J,K
DOUBLE PRECISION TEMP_C,VIRTUAL_T
DOUBLE PRECISION GONV,RONV,SONV
DOUBLE PRECISION FACTOR_G,FACTOR_R,FACTOR_S
DOUBLE PRECISION FACTORB_G,FACTORB_R,FACTORB_S
DOUBLE PRECISION RHOAIR,Z_E
c Constants used to calculate variable intercepts
DOUBLE PRECISION R1,RON,RON2,SON,GON
DOUBLE PRECISION RON_MIN,RON_QR0,RON_DELQR0
DOUBLE PRECISION RON_CONST1R,RON_CONST2R
c Constant intercepts
DOUBLE PRECISION RN0_R,RN0_S,RN0_G
c Other constants
DOUBLE PRECISION RHO_R,RHO_S,RHO_G
DOUBLE PRECISION GAMMA_SEVEN,ALPHA
DOUBLE PRECISION RHOWAT,CELKEL,PI,RD
c Constants used to calculate variable intercepts
R1 = 1.D-15
RON = 8.D6
RON2 = 1.D10
SON = 2.D7
GON = 5.D7
RON_MIN = 8.D6
RON_QR0 = 0.00010D0
RON_DELQR0 = 0.25D0*RON_QR0
RON_CONST1R = (RON2-RON_MIN)*0.5D0
RON_CONST2R = (RON2+RON_MIN)*0.5D0
c Constant intercepts
RN0_R = 8.D6
RN0_S = 2.D7
RN0_G = 4.D6
c Other constants
GAMMA_SEVEN = 720.D0
RHOWAT = 1000.D0
RHO_R = RHOWAT
RHO_S = 100.D0
RHO_G = 400.D0
ALPHA = 0.224D0
CELKEL = 273.15D0
PI = 3.141592653589793D0
RD = 287.04D0
c Force all Q arrays to be 0.0 or greater.
DO K = 1,BTDIM
DO J = 1,SNDIM
DO I = 1,WEDIM
IF (QVP(I,J,K).LT.0.0) THEN
QVP(I,J,K) = 0.0
END IF
IF (QRA(I,J,K).LT.0.0) THEN
QRA(I,J,K) = 0.0
END IF
IF (QSN(I,J,K).LT.0.0) THEN
QSN(I,J,K) = 0.0
END IF
IF (QGR(I,J,K).LT.0.0) THEN
QGR(I,J,K) = 0.0
END IF
END DO
END DO
END DO
c Input pressure is Pa, but we need hPa in calculations
IF (SN0.EQ.0) THEN
DO K = 1,BTDIM
DO J = 1,SNDIM
DO I = 1,WEDIM
IF (TMK(I,J,K).LT.CELKEL) THEN
QSN(I,J,K) = QRA(I,J,K)
QRA(I,J,K) = 0.D0
END IF
END DO
END DO
END DO
END IF
FACTOR_R = GAMMA_SEVEN*1.D18* (1.D0/ (PI*RHO_R))**1.75D0
FACTOR_S = GAMMA_SEVEN*1.D18* (1.D0/ (PI*RHO_S))**1.75D0*
+ (RHO_S/RHOWAT)**2*ALPHA
FACTOR_G = GAMMA_SEVEN*1.D18* (1.D0/ (PI*RHO_G))**1.75D0*
+ (RHO_G/RHOWAT)**2*ALPHA
DO K = 1,BTDIM
DO J = 1,SNDIM
DO I = 1,WEDIM
VIRTUAL_T = TMK(I,J,K)* (0.622D0+QVP(I,J,K))/
+ (0.622D0* (1.D0+QVP(I,J,K)))
RHOAIR = PRS(I,J,K) / (RD*VIRTUAL_T)
c Adjust factor for brightband, where snow or graupel particle
c scatters like liquid water (alpha=1.0) because it is assumed to
c have a liquid skin.
IF (ILIQSKIN.EQ.1 .AND. TMK(I,J,K).GT.CELKEL) THEN
FACTORB_S = FACTOR_S/ALPHA
FACTORB_G = FACTOR_G/ALPHA
ELSE
FACTORB_S = FACTOR_S
FACTORB_G = FACTOR_G
END IF
c Calculate variable intercept parameters
IF (IVARINT.EQ.1) THEN
TEMP_C = DMIN1(-0.001D0,TMK(I,J,K)-CELKEL)
SONV = DMIN1(2.0D8,2.0D6*EXP(-0.12D0*TEMP_C))
GONV = GON
IF (QGR(I,J,K).GT.R1) THEN
GONV = 2.38D0* (PI*RHO_G/
+ (RHOAIR*QGR(I,J,K)))**0.92D0
GONV = MAX(1.D4,MIN(GONV,GON))
END IF
RONV = RON2
IF (QRA(I,J,K).GT.R1) THEN
RONV = RON_CONST1R*TANH((RON_QR0-QRA(I,J,K))/
+ RON_DELQR0) + RON_CONST2R
END IF
ELSE
RONV = RN0_R
SONV = RN0_S
GONV = RN0_G
END IF
c Total equivalent reflectivity factor (z_e, in mm^6 m^-3) is
c the sum of z_e for each hydrometeor species:
Z_E = FACTOR_R* (RHOAIR*QRA(I,J,K))**1.75D0/
+ RONV**.75D0 + FACTORB_S*
+ (RHOAIR*QSN(I,J,K))**1.75D0/SONV**.75D0 +
+ FACTORB_G* (RHOAIR*QGR(I,J,K))**1.75D0/
+ GONV**.75D0
c Adjust small values of Z_e so that dBZ is no lower than -30
Z_E = MAX(Z_E,.001D0)
c Convert to dBZ
DBZ(I,J,K) = 10.D0*LOG10(Z_E)
END DO
END DO
END DO
RETURN
END

511
ncl_reference/wrf_user_latlon_routines.f
Unescape View File

@ -1,511 +0,0 @@
C NCLFORTSTART
SUBROUTINE DLLTOIJ(MAP_PROJ,TRUELAT1,TRUELAT2,STDLON,LAT1,LON1,
+ POLE_LAT,POLE_LON,KNOWNI,KNOWNJ,DX,DY,LATINC,
+ LONINC,LAT,LON,LOC)
DOUBLE PRECISION DELTALON1
DOUBLE PRECISION TL1R
ccc Converts input lat/lon values to the cartesian (i,j) value
ccc for the given projection.
INTEGER MAP_PROJ
DOUBLE PRECISION TRUELAT1,TRUELAT2,STDLON
DOUBLE PRECISION LAT1,LON1,POLE_LAT,POLE_LON,KNOWNI,KNOWNJ
DOUBLE PRECISION DX,DY,LATINC,LONINC,LAT,LON,LOC(2)
C NCLEND
DOUBLE PRECISION CLAIN,DLON,RSW,DELTALON,DELTALAT
DOUBLE PRECISION REFLON,SCALE_TOP,ALA1,ALO1,ALA,ALO,RM,POLEI,POLEJ
c Earth radius divided by dx
DOUBLE PRECISION REBYDX
DOUBLE PRECISION DELTALON1TL1R,CTL1R,ARG,CONE,HEMI
DOUBLE PRECISION I,J
DOUBLE PRECISION LAT1N,LON1N,OLAT,OLON
DOUBLE PRECISION PI,RAD_PER_DEG,DEG_PER_RAD,RE_M
ccc lat1 ! SW latitude (1,1) in degrees (-90->90N)
ccc lon1 ! SW longitude (1,1) in degrees (-180->180E)
ccc dx ! Grid spacing in meters at truelats
ccc dlat ! Lat increment for lat/lon grids
ccc dlon ! Lon increment for lat/lon grids
ccc stdlon ! Longitude parallel to y-axis (-180->180E)
ccc truelat1 ! First true latitude (all projections)
ccc truelat2 ! Second true lat (LC only)
ccc hemi ! 1 for NH, -1 for SH
ccc cone ! Cone factor for LC projections
ccc polei ! Computed i-location of pole point
ccc polej ! Computed j-location of pole point
ccc rsw ! Computed radius to SW corner
ccc knowni ! X-location of known lat/lon
ccc knownj ! Y-location of known lat/lon
ccc RE_M ! Radius of spherical earth, meters
ccc REbydx ! Earth radius divided by dx
PI = 3.141592653589793D0
RAD_PER_DEG = PI/180.D0
DEG_PER_RAD = 180.D0/PI
c Radius of spherical earth, meters
RE_M = 6370000.D0
REBYDX = RE_M/DX
HEMI = 1.0D0
IF (TRUELAT1.LT.0.0D0) THEN
HEMI = -1.0D0
END IF
ccc !MERCATOR
IF (MAP_PROJ.EQ.3) THEN
ccc ! Preliminary variables
CLAIN = COS(RAD_PER_DEG*TRUELAT1)
DLON = DX/ (RE_M*CLAIN)
ccc ! Compute distance from equator to origin, and store in
ccc ! the rsw tag.
RSW = 0.D0
IF (LAT1.NE.0.D0) THEN
RSW = (DLOG(TAN(0.5D0* ((LAT1+90.D0)*RAD_PER_DEG))))/DLON
END IF
DELTALON = LON - LON1
IF (DELTALON.LT.-180.D0) DELTALON = DELTALON + 360.D0
IF (DELTALON.GT.180.D0) DELTALON = DELTALON - 360.D0
I = KNOWNI + (DELTALON/ (DLON*DEG_PER_RAD))
J = KNOWNJ + (DLOG(TAN(0.5D0* ((LAT+90.D0)*RAD_PER_DEG))))/
+ DLON - RSW
ccc !PS
ELSE IF (MAP_PROJ.EQ.2) THEN
REFLON = STDLON + 90.D0
ccc ! Compute numerator term of map scale factor
SCALE_TOP = 1.D0 + HEMI*SIN(TRUELAT1*RAD_PER_DEG)
ccc ! Compute radius to lower-left (SW) corner
ALA1 = LAT1*RAD_PER_DEG
RSW = REBYDX*COS(ALA1)*SCALE_TOP/ (1.D0+HEMI*SIN(ALA1))
ccc ! Find the pole point
ALO1 = (LON1-REFLON)*RAD_PER_DEG
POLEI = KNOWNI - RSW*COS(ALO1)
POLEJ = KNOWNJ - HEMI*RSW*SIN(ALO1)
ccc ! Find radius to desired point
ALA = LAT*RAD_PER_DEG
RM = REBYDX*COS(ALA)*SCALE_TOP/ (1.D0+HEMI*SIN(ALA))
ALO = (LON-REFLON)*RAD_PER_DEG
I = POLEI + RM*COS(ALO)
J = POLEJ + HEMI*RM*SIN(ALO)
ccc !LAMBERT
ELSE IF (MAP_PROJ.EQ.1) THEN
IF (ABS(TRUELAT2).GT.90.D0) THEN
TRUELAT2 = TRUELAT1
END IF
IF (ABS(TRUELAT1-TRUELAT2).GT.0.1D0) THEN
CONE = (DLOG(COS(TRUELAT1*RAD_PER_DEG))-
+ DLOG(COS(TRUELAT2*RAD_PER_DEG)))/
+ (DLOG(TAN((90.D0-ABS(TRUELAT1))*RAD_PER_DEG*
+ 0.5D0))-DLOG(TAN((90.D0-ABS(TRUELAT2))*RAD_PER_DEG*
+ 0.5D0)))
ELSE
CONE = SIN(ABS(TRUELAT1)*RAD_PER_DEG)
END IF
ccc ! Compute longitude differences and ensure we stay
ccc ! out of the forbidden "cut zone"
DELTALON1 = LON1 - STDLON
IF (DELTALON1.GT.+180.D0) DELTALON1 = DELTALON1 - 360.D0
IF (DELTALON1.LT.-180.D0) DELTALON1 = DELTALON1 + 360.D0
ccc ! Convert truelat1 to radian and compute COS for later use
TL1R = TRUELAT1*RAD_PER_DEG
CTL1R = COS(TL1R)
ccc ! Compute the radius to our known lower-left (SW) corner
RSW = REBYDX*CTL1R/CONE* (TAN((90.D0*HEMI-
+ LAT1)*RAD_PER_DEG/2.D0)/TAN((90.D0*HEMI-
+ TRUELAT1)*RAD_PER_DEG/2.D0))**CONE
ccc ! Find pole point
ARG = CONE* (DELTALON1*RAD_PER_DEG)
POLEI = HEMI*KNOWNI - HEMI*RSW*SIN(ARG)
POLEJ = HEMI*KNOWNJ + RSW*COS(ARG)
ccc ! Compute deltalon between known longitude and standard
ccc ! lon and ensure it is not in the cut zone
DELTALON = LON - STDLON
IF (DELTALON.GT.+180.D0) DELTALON = DELTALON - 360.D0
IF (DELTALON.LT.-180.D0) DELTALON = DELTALON + 360.D0
ccc ! Radius to desired point
RM = REBYDX*CTL1R/CONE* (TAN((90.D0*HEMI-
+ LAT)*RAD_PER_DEG/2.D0)/TAN((90.D0*HEMI-
+ TRUELAT1)*RAD_PER_DEG/2.D0))**CONE
ARG = CONE* (DELTALON*RAD_PER_DEG)
I = POLEI + HEMI*RM*SIN(ARG)
J = POLEJ - RM*COS(ARG)
ccc ! Finally, if we are in the southern hemisphere, flip the
ccc ! i/j values to a coordinate system where (1,1) is the SW
ccc ! corner (what we assume) which is different than the
ccc ! original NCEP algorithms which used the NE corner as
ccc ! the origin in the southern hemisphere (left-hand vs.
ccc ! right-hand coordinate?)
I = HEMI*I
J = HEMI*J
ccc !lat-lon
ELSE IF (MAP_PROJ.EQ.6) THEN
IF (POLE_LAT.NE.90.D0) THEN
CALL ROTATECOORDS(LAT,LON,OLAT,OLON,POLE_LAT,POLE_LON,
+ STDLON,-1)
LAT = OLAT
LON = OLON + STDLON
END IF
c ! make sure center lat/lon is good
IF (POLE_LAT.NE.90.D0) THEN
CALL ROTATECOORDS(LAT1,LON1,OLAT,OLON,POLE_LAT,POLE_LON,
+ STDLON,-1)
LAT1N = OLAT
LON1N = OLON + STDLON
DELTALAT = LAT - LAT1N
DELTALON = LON - LON1N
ELSE
DELTALAT = LAT - LAT1
DELTALON = LON - LON1
END IF
c ! Compute i/j
I = DELTALON/LONINC
J = DELTALAT/LATINC
I = I + KNOWNI
J = J + KNOWNJ
ELSE
PRINT *,'ERROR: Do not know map projection ',MAP_PROJ
END IF
LOC(1) = J
LOC(2) = I
RETURN
END
C NCLFORTSTART
SUBROUTINE DIJTOLL(MAP_PROJ,TRUELAT1,TRUELAT2,STDLON,LAT1,LON1,
+ POLE_LAT,POLE_LON,KNOWNI,KNOWNJ,DX,DY,LATINC,
+ LONINC,AI,AJ,LOC)
DOUBLE PRECISION GI2
DOUBLE PRECISION ARCCOS
DOUBLE PRECISION DELTALON1
DOUBLE PRECISION TL1R
ccc ! Converts input lat/lon values to the cartesian (i,j) value
ccc ! for the given projection.
INTEGER MAP_PROJ
DOUBLE PRECISION TRUELAT1,TRUELAT2,STDLON
DOUBLE PRECISION LAT1,LON1,POLE_LAT,POLE_LON,KNOWNI,KNOWNJ
DOUBLE PRECISION DX,DY,LATINC,LONINC,AI,AJ,LOC(2)
C NCLEND
DOUBLE PRECISION CLAIN,DLON,RSW,DELTALON,DELTALAT
DOUBLE PRECISION REFLON,SCALE_TOP,ALA1,ALO1,ALA,ALO,RM,POLEI,POLEJ
c Earth radius divided by dx
DOUBLE PRECISION REBYDX
DOUBLE PRECISION DELTALON1TL1R,CTL1R,ARG,CONE,HEMI
DOUBLE PRECISION PI,RAD_PER_DEG,DEG_PER_RAD,RE_M
DOUBLE PRECISION INEW,JNEW,R,R2
DOUBLE PRECISION CHI,CHI1,CHI2
DOUBLE PRECISION XX,YY,LAT,LON
DOUBLE PRECISION RLAT,RLON,OLAT,OLON,LAT1N,LON1N
DOUBLE PRECISION PHI_NP,LAM_NP,LAM_0,DLAM
DOUBLE PRECISION SINPHI,COSPHI,COSLAM,SINLAM
ccc lat1 ! SW latitude (1,1) in degrees (-90->90N)
ccc lon1 ! SW longitude (1,1) in degrees (-180->180E)
ccc dx ! Grid spacing in meters at truelats
ccc dlat ! Lat increment for lat/lon grids
ccc dlon ! Lon increment for lat/lon grids
ccc stdlon ! Longitude parallel to y-axis (-180->180E)
ccc truelat1 ! First true latitude (all projections)
ccc truelat2 ! Second true lat (LC only)
ccc hemi ! 1 for NH, -1 for SH
ccc cone ! Cone factor for LC projections
ccc polei ! Computed i-location of pole point
ccc polej ! Computed j-location of pole point
ccc rsw ! Computed radius to SW corner
ccc knowni ! X-location of known lat/lon
ccc knownj ! Y-location of known lat/lon
ccc RE_M ! Radius of spherical earth, meters
ccc REbydx ! Earth radius divided by dx
PI = 3.141592653589793D0
RAD_PER_DEG = PI/180.D0
DEG_PER_RAD = 180.D0/PI
c Radius of spherical earth, meters
RE_M = 6370000.D0
REBYDX = RE_M/DX
HEMI = 1.0D0
IF (TRUELAT1.LT.0.0D0) THEN
HEMI = -1.0D0
END IF
ccc !MERCATOR
IF (MAP_PROJ.EQ.3) THEN
ccc ! Preliminary variables
CLAIN = COS(RAD_PER_DEG*TRUELAT1)
DLON = DX/ (RE_M*CLAIN)
ccc ! Compute distance from equator to origin, and store in
ccc ! the rsw tag.
RSW = 0.D0
IF (LAT1.NE.0.D0) THEN
RSW = (DLOG(TAN(0.5D0* ((LAT1+90.D0)*RAD_PER_DEG))))/DLON
END IF
LAT = 2.0D0*ATAN(EXP(DLON* (RSW+AJ-KNOWNJ)))*DEG_PER_RAD -
+ 90.D0
LON = (AI-KNOWNI)*DLON*DEG_PER_RAD + LON1
IF (LON.GT.180.D0) LON = LON - 360.D0
IF (LON.LT.-180.D0) LON = LON + 360.D0
ccc !PS
ELSE IF (MAP_PROJ.EQ.2) THEN
ccc ! Compute the reference longitude by rotating 90 degrees to
ccc ! the east to find the longitude line parallel to the
ccc ! positive x-axis.
REFLON = STDLON + 90.D0
ccc ! Compute numerator term of map scale factor
SCALE_TOP = 1.D0 + HEMI*SIN(TRUELAT1*RAD_PER_DEG)
ccc ! Compute radius to known point
ALA1 = LAT1*RAD_PER_DEG
RSW = REBYDX*COS(ALA1)*SCALE_TOP/ (1.D0+HEMI*SIN(ALA1))
ccc ! Find the pole point
ALO1 = (LON1-REFLON)*RAD_PER_DEG
POLEI = KNOWNI - RSW*COS(ALO1)
POLEJ = KNOWNJ - HEMI*RSW*SIN(ALO1)
ccc ! Compute radius to point of interest
XX = AI - POLEI
YY = (AJ-POLEJ)*HEMI
R2 = XX**2 + YY**2
ccc ! Now the magic code
IF (R2.EQ.0.D0) THEN
LAT = HEMI*90.D0
LON = REFLON
ELSE
GI2 = (REBYDX*SCALE_TOP)**2.D0
LAT = DEG_PER_RAD*HEMI*ASIN((GI2-R2)/ (GI2+R2))
ARCCOS = ACOS(XX/SQRT(R2))
IF (YY.GT.0) THEN
LON = REFLON + DEG_PER_RAD*ARCCOS
ELSE
LON = REFLON - DEG_PER_RAD*ARCCOS
END IF
END IF
ccc ! Convert to a -180 -> 180 East convention
IF (LON.GT.180.D0) LON = LON - 360.D0
IF (LON.LT.-180.D0) LON = LON + 360.D0
ccc !LAMBERT
ELSE IF (MAP_PROJ.EQ.1) THEN
IF (ABS(TRUELAT2).GT.90.D0) THEN
TRUELAT2 = TRUELAT1
END IF
IF (ABS(TRUELAT1-TRUELAT2).GT.0.1D0) THEN
CONE = (DLOG(COS(TRUELAT1*RAD_PER_DEG))-
+ DLOG(COS(TRUELAT2*RAD_PER_DEG)))/
+ (DLOG(TAN((90.D0-ABS(TRUELAT1))*RAD_PER_DEG*
+ 0.5D0))-DLOG(TAN((90.D0-ABS(TRUELAT2))*RAD_PER_DEG*
+ 0.5D0)))
ELSE
CONE = SIN(ABS(TRUELAT1)*RAD_PER_DEG)
END IF
ccc ! Compute longitude differences and ensure we stay out of the
ccc ! forbidden "cut zone"
DELTALON1 = LON1 - STDLON
IF (DELTALON1.GT.+180.D0) DELTALON1 = DELTALON1 - 360.D0
IF (DELTALON1.LT.-180.D0) DELTALON1 = DELTALON1 + 360.D0
ccc ! Convert truelat1 to radian and compute COS for later use
TL1R = TRUELAT1*RAD_PER_DEG
CTL1R = COS(TL1R)
ccc ! Compute the radius to our known point
RSW = REBYDX*CTL1R/CONE* (TAN((90.D0*HEMI-
+ LAT1)*RAD_PER_DEG/2.D0)/TAN((90.D0*HEMI-
+ TRUELAT1)*RAD_PER_DEG/2.D0))**CONE
ccc ! Find pole point
ALO1 = CONE* (DELTALON1*RAD_PER_DEG)
POLEI = HEMI*KNOWNI - HEMI*RSW*SIN(ALO1)
POLEJ = HEMI*KNOWNJ + RSW*COS(ALO1)
CHI1 = (90.D0-HEMI*TRUELAT1)*RAD_PER_DEG
CHI2 = (90.D0-HEMI*TRUELAT2)*RAD_PER_DEG
ccc ! See if we are in the southern hemispere and flip the
ccc ! indices if we are.
INEW = HEMI*AI
JNEW = HEMI*AJ
ccc ! Compute radius**2 to i/j location
REFLON = STDLON + 90.D0
XX = INEW - POLEI
YY = POLEJ - JNEW
R2 = (XX*XX+YY*YY)
R = SQRT(R2)/REBYDX
ccc ! Convert to lat/lon
IF (R2.EQ.0.D0) THEN
LAT = HEMI*90.D0
LON = STDLON
ELSE
LON = STDLON + DEG_PER_RAD*ATAN2(HEMI*XX,YY)/CONE
LON = DMOD(LON+360.D0,360.D0)
IF (CHI1.EQ.CHI2) THEN
CHI = 2.0D0*ATAN((R/TAN(CHI1))** (1.D0/CONE)*
+ TAN(CHI1*0.5D0))
ELSE
CHI = 2.0D0*ATAN((R*CONE/SIN(CHI1))** (1.D0/CONE)*
+ TAN(CHI1*0.5D0))
END IF
LAT = (90.0D0-CHI*DEG_PER_RAD)*HEMI
END IF
IF (LON.GT.+180.D0) LON = LON - 360.D0
IF (LON.LT.-180.D0) LON = LON + 360.D0
ccc !lat-lon
ELSE IF (MAP_PROJ.EQ.6) THEN
INEW = AI - KNOWNI
JNEW = AJ - KNOWNJ
IF (INEW.LT.0.D0) INEW = INEW + 360.D0/LONINC
IF (INEW.GE.360.D0/DX) INEW = INEW - 360.D0/LONINC
c
ccc ! Compute deltalat and deltalon
DELTALAT = JNEW*LATINC
DELTALON = INEW*LONINC
IF (POLE_LAT.NE.90.D0) THEN
CALL ROTATECOORDS(LAT1,LON1,OLAT,OLON,POLE_LAT,POLE_LON,
+ STDLON,-1)
LAT1N = OLAT
LON1N = OLON + STDLON
LAT = DELTALAT + LAT1N
LON = DELTALON + LON1N
ELSE
LAT = DELTALAT + LAT1
LON = DELTALON + LON1
END IF
IF (POLE_LAT.NE.90.D0) THEN
LON = LON - STDLON
CALL ROTATECOORDS(LAT,LON,OLAT,OLON,POLE_LAT,POLE_LON,
+ STDLON,1)
LAT = OLAT
LON = OLON
END IF
IF (LON.LT.-180.D0) LON = LON + 360.D0
IF (LON.GT.180.D0) LON = LON - 360.D0
ELSE
PRINT *,'ERROR: Do not know map projection ',MAP_PROJ
END IF
LOC(1) = LAT
LOC(2) = LON
RETURN
END
C NCLFORTSTART
SUBROUTINE ROTATECOORDS(ILAT,ILON,OLAT,OLON,LAT_NP,LON_NP,LON_0,
+ DIRECTION)
DOUBLE PRECISION ILAT,ILON
DOUBLE PRECISION OLAT,OLON
DOUBLE PRECISION LAT_NP,LON_NP,LON_0
INTEGER DIRECTION
C NCLEND
c ! >=0, default : computational -> geographical
c ! < 0 : geographical -> computational
DOUBLE PRECISION RLAT,RLON
DOUBLE PRECISION PHI_NP,LAM_NP,LAM_0,DLAM
DOUBLE PRECISION SINPHI,COSPHI,COSLAM,SINLAM
DOUBLE PRECISION PI,RAD_PER_DEG,DEG_PER_RAD
PI = 3.141592653589793D0
RAD_PER_DEG = PI/180.D0
DEG_PER_RAD = 180.D0/PI
c ! Convert all angles to radians
PHI_NP = LAT_NP*RAD_PER_DEG
LAM_NP = LON_NP*RAD_PER_DEG
LAM_0 = LON_0*RAD_PER_DEG
RLAT = ILAT*RAD_PER_DEG
RLON = ILON*RAD_PER_DEG
IF (DIRECTION.LT.0) THEN
c ! The equations are exactly the same except for one
c ! small difference with respect to longitude ...
DLAM = PI - LAM_0
ELSE
DLAM = LAM_NP
END IF
SINPHI = COS(PHI_NP)*COS(RLAT)*COS(RLON-DLAM) +
+ SIN(PHI_NP)*SIN(RLAT)
COSPHI = SQRT(1.D0-SINPHI*SINPHI)
COSLAM = SIN(PHI_NP)*COS(RLAT)*COS(RLON-DLAM) -
+ COS(PHI_NP)*SIN(RLAT)
SINLAM = COS(RLAT)*SIN(RLON-DLAM)
IF (COSPHI.NE.0.D0) THEN
COSLAM = COSLAM/COSPHI
SINLAM = SINLAM/COSPHI
END IF
OLAT = DEG_PER_RAD*ASIN(SINPHI)
OLON = DEG_PER_RAD* (ATAN2(SINLAM,COSLAM)-DLAM-LAM_0+LAM_NP)
END

404
ncl_reference/wrf_vinterp.f
Unescape View File

@ -1,404 +0,0 @@
CThe subroutines in this file were taken directly from RIP code written
C by Dr. Mark Stoelinga. They were modified by Sherrie
C Fredrick(NCAR/MMM) to work with NCL February 2015.
C NCLFORTSTART
subroutine wrf_monotonic(out,in,lvprs,cor,idir,delta,
& ew,ns,nz,icorsw)
implicit none
integer idir,ew,ns,nz,icorsw
double precision delta
double precision in(ew,ns,nz),out(ew,ns,nz)
double precision lvprs(ew,ns,nz),cor(ew,ns)
C NCLEND
integer i,j,k,ripk,k300
do j=1,ns
do i=1,ew
if (icorsw.eq.1.and.cor(i,j).lt.0.) then
do k=1,nz
in(i,j,k)=-in(i,j,k)
enddo
endif
c
c First find k index that is at or below (height-wise) the 300 hPa
c level.
c
do k = 1,nz
ripk = nz-k+1
if (lvprs(i,j,k) .le. 300.d0) then
k300=k
goto 40
endif
enddo
c
40 continue
do k = k300, 1,-1
if (idir.eq.1) then
out(i,j,k)=min(in(i,j,k),in(i,j,k+1)+delta)
elseif (idir.eq.-1) then
out(i,j,k)=max(in(i,j,k),in(i,j,k+1)-delta)
endif
enddo
do k = k300+1, nz
if (idir.eq.1) then
out(i,j,k)=max(in(i,j,k),in(i,j,k-1)-delta)
elseif (idir.eq.-1) then
out(i,j,k)=min(in(i,j,k),in(i,j,k-1)+delta)
endif
enddo
end do
end do
return
end
c--------------------------------------------------------------------
C NCLFORTSTART
FUNCTION wrf_intrp_value (wvalp0,wvalp1,vlev,vcp0,vcp1,icase)
implicit none
integer icase
double precision wvalp0,wvalp1,vlev,vcp0,vcp1
C NCLEND
double precision valp0,valp1,rvalue,rgas,ussalr,sclht
double precision wrf_intrp_value,chkdiff
rgas = 287.04d0 !J/K/kg
ussalr = 0.0065d0 ! deg C per m
sclht = rgas*256.d0/9.81d0
valp0 = wvalp0
valp1 = wvalp1
if ( icase .eq. 2) then !GHT
valp0=exp(-wvalp0/sclht)
valp1=exp(-wvalp1/sclht)
end if
chkdiff = vcp1 - vcp0
if(chkdiff .eq. 0) then
print *,"bad difference in vcp's"
stop
end if
rvalue = (vlev-vcp0)*(valp1-valp0)/(vcp1-vcp0)+valp0
if (icase .eq. 2) then !GHT
wrf_intrp_value = -sclht*log(rvalue)
else
wrf_intrp_value = rvalue
endif
return
end
c------------------------------------------------------------
C NOTES:
c vcarray is the array holding the values for the vertical
c coordinate.
c It will always come in with the dimensions of
c the staggered U and V grid.
C NCLFORTSTART
subroutine wrf_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,
& sfp,smsfp,vcarray,interp_levels,numlevels,
& icase,ew,ns,nz,extrap,vcor,logp,rmsg)
implicit none
integer ew,ns,nz,icase,extrap
integer vcor,numlevels,logp
double precision datain(ew,ns,nz),pres(ew,ns,nz),tk(ew,ns,nz)
double precision ght(ew,ns,nz)
double precision terrain(ew,ns),sfp(ew,ns),smsfp(ew,ns)
double precision dataout(ew,ns,numlevels),qvp(ew,ns,nz)
double precision vcarray(ew,ns,nz)
double precision interp_levels(numlevels),rmsg
C NCLEND
integer njx,niy,nreqlvs,ripk
integer i,j,k,itriv,kupper
integer ifound,miy,mjx,isign
double precision rlevel,vlev,diff
double precision tempout(ew,ns),tmpvlev
double precision vcp1,vcp0,valp0,valp1
double precision rgas,rgasmd,sclht,ussalr,cvc,eps
double precision qvlhsl,ttlhsl,vclhsl,vctophsl
double precision wrf_intrp_value
double precision plhsl,zlhsl,ezlhsl,tlhsl,psurf,pratio,tlev
double precision ezsurf,psurfsm,zsurf,qvapor,vt
double precision rconst,expon,exponi
double precision ezlev,plev,zlev,ptarget,dpmin,dp
double precision pbot,zbot,tbotextrap,e
double precision tlclc1,tlclc2,tlclc3,tlclc4
double precision thtecon1,thtecon2,thtecon3
double precision tlcl,gamma,cp,cpmd,gammamd,gammam
character cvcord*1
rgas = 287.04d0 !J/K/kg
rgasmd = .608d0
ussalr = .0065d0 ! deg C per m
sclht = rgas*256.d0/9.81d0
eps = 0.622d0
rconst = -9.81d0/(rgas * ussalr)
expon = rgas*ussalr/9.81d0
exponi = 1./expon
tlclc1 = 2840.d0
tlclc2 = 3.5d0
tlclc3 = 4.805d0
tlclc4 = 55.d0
thtecon1 = 3376.d0 ! K
thtecon2 = 2.54d0
thtecon3 = 0.81d0
cp = 1004.d0
cpmd = 0.887d0
gamma = rgas/cp
gammamd = rgasmd-cpmd
if(vcor .eq. 1) then
cvcord = 'p'
else if((vcor .eq. 2) .or. (vcor .eq. 3)) then
cvcord = 'z'
else if((vcor .eq. 4) .or. (vcor .eq. 5)) then
cvcord = 't'
end if
miy = ns
mjx = ew
njx = ew
niy = ns
do j = 1,mjx
do i = 1,miy
tempout(j,i) = rmsg
end do
end do
do nreqlvs = 1,numlevels
if(cvcord .eq. 'z') then
!Convert rlevel to meters from km
rlevel = interp_levels(nreqlvs) * 1000.d0
vlev = exp(-rlevel/sclht)
else if(cvcord .eq. 'p') then
vlev = interp_levels(nreqlvs)
else if(cvcord .eq. 't') then
vlev = interp_levels(nreqlvs)
end if
do j=1,mjx
do i=1,miy
cGet the interpolated value that is within the model domain
ifound = 0
do k = 1,nz-1
ripk = nz-k+1
vcp1 = vcarray(j,i,ripk-1)
vcp0 = vcarray(j,i,ripk)
valp0 = datain(j,i,ripk)
valp1 = datain(j,i,ripk-1)
if ((vlev.ge.vcp0.and.vlev.le.vcp1) .or.
& (vlev.le.vcp0.and.vlev.ge.vcp1)) then
c print *,i,j,valp0,valp1
if((valp0 .eq. rmsg).or.(valp1 .eq. rmsg)) then
tempout(j,i) = rmsg
ifound=1
else
if(logp .eq. 1) then
vcp1 = log(vcp1)
vcp0 = log(vcp0)
if(vlev .eq. 0.0d0) then
print *,"Pressure value = 0"
print *,"Unable to take log of 0"
stop
end if
tmpvlev = log(vlev)
else
tmpvlev = vlev
end if
tempout(j,i) = wrf_intrp_value(valp0,valp1,
& tmpvlev,vcp0,vcp1,icase)
c print *,"one ",i,j,tempout(j,i)
ifound=1
end if
goto 115
end if
end do !end for the k loop
115 continue
if (ifound.eq.1) then !Grid point is in the model domain
goto 333
end if
cIf the user has requested no extrapolatin then just assign
call values above or below the model level to rmsg.
if(extrap .eq. 0) then
tempout(j,i) = rmsg
goto 333
end if
c The grid point is either above or below the model domain
c
c First we will check to see if the grid point is above the
c model domain.
vclhsl = vcarray(j,i,1) !lowest model level
vctophsl = vcarray(j,i,nz)!highest model level
diff = vctophsl-vclhsl
isign = nint(diff/abs(diff))
C
if(isign*vlev.ge.isign*vctophsl) then
C Assign the highest model level to the out array
tempout(j,i)=datain(j,i,nz)
C print *,"at warn",j,i,tempout(j,i)
goto 333
endif
c
c Only remaining possibility is that the specified level is below
c lowest model level. If lowest model level value is missing,
c set interpolated value to missing.
c
if (datain(i,j,1) .eq. rmsg) then
tempout(j,i) = rmsg
goto 333
endif
c
c If the field comming in is not a pressure,temperature or height
C field we can set the output array to the value at the lowest
c model level.
c
tempout(j,i) = datain(j,i,1)
c
c For the special cases of pressure on height levels or height on
c pressure levels, or temperature-related variables on pressure or
c height levels, perform a special extrapolation based on
c US Standard Atmosphere. Here we calcualate the surface pressure
c with the altimeter equation. This is how RIP calculates the
c surface pressure.
c
if (icase.gt.0) then
plhsl = pres(j,i,1) * 0.01d0 !pressure at lowest model level
zlhsl = ght(j,i,1) !grid point height a lowest model level
ezlhsl = exp(-zlhsl/sclht)
tlhsl = tk(j,i,1) !temperature in K at lowest model level
zsurf = terrain(j,i)
qvapor = max((qvp(j,i,1)*.001d0),1.e-15)
c virtual temperature
c vt = tlhsl * (eps + qvapor)/(eps*(1.0 + qvapor))
c psurf = plhsl * (vt/(vt+ussalr * (zlhsl-zsurf)))**rconst
psurf = sfp(j,i)
psurfsm = smsfp(j,i)
ezsurf = exp(-zsurf/sclht)
cThe if for checking above ground
if ((cvcord.eq.'z'.and.vlev.lt.ezsurf).or.
& (cvcord.eq.'p'.and.vlev.lt.psurf)) then
c
c We are below the lowest data level but above the ground.
c Use linear interpolation (linear in prs and exp-height).
c
if (cvcord.eq.'p') then
plev=vlev
ezlev=((plev-plhsl)*ezsurf+(psurf-plev)*ezlhsl)/
& (psurf-plhsl)
zlev=-sclht*log(ezlev)
if (icase .eq. 2) then
tempout(j,i)=zlev
goto 333
endif
elseif (cvcord.eq.'z') then
ezlev=vlev
zlev=-sclht*log(ezlev)
plev=((ezlev-ezlhsl)*psurf+(ezsurf-ezlev)*plhsl)/
& (ezsurf-ezlhsl)
if (icase .eq. 1) then
tempout(j,i)=plev
goto 333
endif
endif
else !else for checking above ground
ptarget=psurfsm-150.d0
dpmin=1.e4
do k=1,nz
ripk = nz-k+1
dp=abs((pres(j,i,ripk) * 0.01d0)-ptarget)
if (dp.gt.dpmin) goto 334
dpmin=min(dpmin,dp)
enddo
334 kupper=k-1
ripk = nz - kupper + 1
pbot = max(plhsl,psurf)
zbot = min(zlhsl,zsurf)
pratio = pbot/(pres(j,i,ripk) * 0.01d0)
tbotextrap = tk(j,i,ripk)*(pratio)**expon
c virtual temperature
vt = tbotextrap * (eps + qvapor)/(eps*(1.0d0+qvapor))
if (cvcord.eq.'p') then
plev=vlev
zlev=zbot+vt/ussalr*(1.-(vlev/pbot)**expon)
if(icase .eq. 2) then
tempout(j,i)=zlev
goto 333
endif
elseif (cvcord.eq.'z') then
zlev=-sclht*log(vlev)
plev=pbot*(1.+ussalr/vt*(zbot-zlev))**exponi
if (icase .eq. 1) then
tempout(j,i)=plev
goto 333
endif
endif
end if !end if for checking above ground
end if !for icase gt 0
if(icase .gt. 2) then !extrapolation for temperature
tlev=tlhsl+(zlhsl-zlev)*ussalr
qvapor = max(qvp(j,i,1),1.e-15)
gammam = gamma*(1.+gammamd*qvapor)
if(icase .eq. 3) then
tempout(j,i) = tlev - 273.16d0
else if(icase .eq. 4) then
tempout(j,i) = tlev
C Potential temperature - theta
else if (icase. eq. 5) then
tempout(j,i)=tlev*(1000.d0/plev)**gammam
C extraolation for equivalent potential temperature
else if (icase .eq. 6) then
e = qvapor*plev/(eps+qvapor)
tlcl = tlclc1/(log(tlev**tlclc2/e)-tlclc3)+tlclc4
tempout(j,i)=tlev*(1000.d0/plev)**(gammam)*
& exp((thtecon1/tlcl-thtecon2)*qvapor*
& (1.+thtecon3*qvapor))
end if
end if
333 continue
end do
end do
! print *,"----done----",interp_levels(nreqlvs)
do i = 1,njx
do j = 1,niy
dataout(i,j,nreqlvs) = tempout(i,j)
end do
end do
end do !end for the nreqlvs
return
end !wrf_vinterp

1147
ncl_reference/wrf_vinterpW.c
View File

File diff suppressed because it is too large Load Diff
Write Preview
Loading…
Cancel
Save