SUBROUTINE f_interpz3d(data3d,zdata,desiredloc,missingval,out2d,nx,ny,nz)

IMPLICIT NONE

INTEGER, INTENT(IN) :: nx,ny,nz

REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: data3d

REAL(KIND=8), DIMENSION(nx,ny), INTENT(OUT) :: out2d

REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: zdata

REAL(KIND=8), INTENT(IN) :: desiredloc

REAL(KIND=8), INTENT(IN) :: missingval

INTEGER :: i,j,kp,ip,im

LOGICAL :: dointerp

REAL(KIND=8) :: height,w1,w2

height = desiredloc

! does vertical coordinate increase or decrease with increasing k?

! set offset appropriately

ip = 0

im = 1

IF (zdata(1,1,1).GT.zdata(1,1,nz)) THEN

ip = 1

im = 0

END IF

DO i = 1,nx

DO j = 1,ny

! Initialize to missing. Was initially hard-coded to -999999.

out2d(i,j) = missingval

dointerp = .FALSE.

kp = nz

DO WHILE ((.NOT. dointerp) .AND. (kp >= 2))

IF (((zdata(i,j,kp-im) < height) .AND. (zdata(i,j,kp-ip) > height))) THEN

w2 = (height-zdata(i,j,kp-im))/(zdata(i,j,kp-ip)-zdata(i,j,kp-im))

w1 = 1.D0 - w2

out2d(i,j) = w1*data3d(i,j,kp-im) + w2*data3d(i,j,kp-ip)

dointerp = .TRUE.

END IF

kp = kp - 1

END DO

END DO

END DO

RETURN

END SUBROUTINE f_interpz3d

SUBROUTINE f_interp2dxy(v3d,xy,v2d,nx,ny,nz,nxy)

IMPLICIT NONE

INTEGER,INTENT(IN) :: nx,ny,nz,nxy

REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(IN) :: v3d

REAL(KIND=8),DIMENSION(nxy,nz),INTENT(OUT) :: v2d

REAL(KIND=8),DIMENSION(2,nxy),INTENT(IN) :: xy

INTEGER :: i,j,k,ij

REAL(KIND=8) :: 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 SUBROUTINE f_interp2dxy

SUBROUTINE f_interp1d(v_in,z_in,z_out,vmsg,v_out,nz_in,nz_out)

IMPLICIT NONE

INTEGER,INTENT(IN) :: nz_in, nz_out

REAL(KIND=8),DIMENSION(nz_in),INTENT(IN) :: v_in,z_in

REAL(KIND=8),DIMENSION(nz_out),INTENT(IN) :: z_out

REAL(KIND=8),DIMENSION(nz_out),INTENT(OUT) :: v_out

REAL(KIND=8),INTENT(IN) :: vmsg

INTEGER :: kp,k,im,ip

LOGICAL :: interp

REAL(KIND=8) :: height,w1,w2

! does vertical coordinate increase of decrease with increasing k?

! 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 SUBROUTINE f_interp1d

! This routine assumes

! index order is (i,j,k)

! wrf staggering

!

! units: pressure (Pa), temperature(K), height (m), mixing ratio

! (kg kg{-1}) availability of 3d p, t, and qv; 2d terrain; 1d

! half-level zeta string

! output units of SLP are Pa, but you should divide that by 100 for the

! weather weenies.

! virtual effects are included

!

SUBROUTINE f_computeslp(z,t,p,q,t_sea_level,t_surf,level,throw_exception,&

sea_level_pressure,nx,ny,nz)

IMPLICIT NONE

EXTERNAL throw_exception

! Estimate sea level pressure.

INTEGER, INTENT(IN) :: nx,ny,nz

REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(IN) :: z

REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(IN) :: t,p,q

! The output is the 2d sea level pressure.

REAL(KIND=8),DIMENSION(nx,ny),INTENT(OUT) :: sea_level_pressure

INTEGER,DIMENSION(nx,ny), INTENT(INOUT) :: level

REAL(KIND=8),DIMENSION(nx,ny),INTENT(INOUT) :: t_surf,t_sea_level

! Some required physical constants:

REAL(KIND=8), PARAMETER :: R=287.04D0, G=9.81D0, GAMMA=0.0065D0

! Specific constants for assumptions made in this routine:

REAL(KIND=8), PARAMETER :: TC=273.16D0+17.5D0, PCONST=10000

LOGICAL, PARAMETER :: ridiculous_mm5_test=.TRUE.

! PARAMETER (ridiculous_mm5_test = .FALSE.)

! Local variables:

INTEGER :: i,j,k

INTEGER :: klo,khi

REAL(KIND=8) :: plo,phi,tlo,thi,zlo,zhi

REAL(KIND=8) :: p_at_pconst,t_at_pconst,z_at_pconst

REAL(KIND=8) :: z_half_lowest

LOGICAL :: l1,l2,l3,found

! Find least zeta level that is PCONST Pa above the surface. We

! later use this level to extrapolate a surface pressure and

! temperature, which is supposed to reduce the effect of the diurnal

! 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 <= nz))

IF (p(i,j,k) < p(i,j,1)-PCONST) THEN

level(i,j) = k

found = .TRUE.

END IF

k = k + 1

END DO

IF (level(i,j) == -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.'

CALL throw_exception('Error in finding 100 hPa up')

END IF

END DO

END DO

! Get temperature PCONST Pa above surface. Use this to extrapolate

! 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 == khi) THEN

!PRINT '(A)','Trapping levels are weird.'

!PRINT '(A,I3,A,I3,A)','klo = ',klo,', khi = ',khi,': and they should not be equal.'

CALL throw_exception('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))

! zlo = zetahalf(klo)/ztop*(ztop-terrain(i,j))+terrain(i,j)

! 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

! If we follow a traditional computation, there is a correction to the

! sea level temperature if both the surface and sea level

! temperatures are *too* hot.

IF (ridiculous_mm5_test) THEN

DO J = 1,ny

DO I = 1,nx

l1 = t_sea_level(i,j) < TC

l2 = t_surf(i,j) <= 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

! The grand finale: ta da!

DO J = 1,ny

DO I = 1,nx

! z_half_lowest=zetahalf(1)/ztop*(ztop-terrain(i,j))+terrain(i,j)

z_half_lowest = z(i,j,1)

! Convert to hPa in this step, by multiplying by 0.01. The original

! Fortran routine didn't do this, but the NCL script that called it

! 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

! PRINT *,'sea pres input at weird location i=20,j=1,k=1'

! PRINT *,'t=',t(20,1,1),t(20,2,1),t(20,3,1)

! PRINT *,'z=',z(20,1,1),z(20,2,1),z(20,3,1)

! PRINT *,'p=',p(20,1,1),p(20,2,1),p(20,3,1)

! PRINT *,'slp=',sea_level_pressure(20,1),

! * sea_level_pressure(20,2),sea_level_pressure(20,3)

RETURN

END SUBROUTINE f_computeslp

! Temperature from potential temperature in kelvin.

SUBROUTINE f_computetk(pressure,theta,tk,nx)

IMPLICIT NONE

INTEGER, INTENT(IN) :: nx

REAL(KIND=8) :: pi

REAL(KIND=8), DIMENSION(nx), INTENT(IN) :: pressure

REAL(KIND=8), DIMENSION(nx), INTENT(IN) :: theta

REAL(KIND=8), DIMENSION(nx), INTENT(OUT) :: tk

INTEGER :: i

REAL(KIND=8), 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

RETURN

END SUBROUTINE f_computetk

! Dewpoint. Note: 1D array arguments.

SUBROUTINE f_computetd(pressure,qv_in,td,nx)

IMPLICIT NONE

INTEGER, INTENT(IN) :: nx

REAL(KIND=8),DIMENSION(nx),INTENT(IN) :: pressure

REAL(KIND=8),DIMENSION(nx),INTENT(IN) :: qv_in

REAL(KIND=8),DIMENSION(nx),INTENT(OUT) :: td

REAL(KIND=8) :: qv,tdc

INTEGER :: i

DO i = 1,nx

qv = DMAX1(qv_in(i),0.D0)

! vapor pressure

tdc = qv*pressure(i)/ (.622D0 + qv)

! 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 SUBROUTINE f_computetd

! Relative Humidity. Note: 1D array arguments

SUBROUTINE f_computerh(qv,p,t,rh,nx)

IMPLICIT NONE

INTEGER, INTENT(IN) :: nx

REAL(KIND=8),DIMENSION(nx),INTENT(IN) :: qv,p,t

REAL(KIND=8),DIMENSION(nx),INTENT(OUT) :: rh

REAL(KIND=8), PARAMETER :: SVP1=0.6112D0,SVP2=17.67D0,SVP3=29.65D0,SVPT0=273.15D0

INTEGER :: i

REAL(KIND=8) :: qvs,es,pressure,temperature

REAL(KIND=8), PARAMETER :: R_D=287.D0,R_V=461.6D0,EP_2=R_D/R_V

REAL(KIND=8), PARAMETER :: EP_3=0.622D0

DO i = 1,nx

pressure = p(i)

temperature = t(i)

! es = 1000.*svp1*

es = 10.D0*SVP1*EXP(SVP2* (temperature-SVPT0)/(temperature-SVP3))

! qvs = ep_2*es/(pressure-es)

qvs = EP_3*es/ (0.01D0*pressure- (1.D0-EP_3)*es)

! rh = 100*amax1(1., qv(i)/qvs)

! rh(i) = 100.*qv(i)/qvs

rh(i) = 100.D0*DMAX1(DMIN1(qv(i)/qvs,1.0D0),0.0D0)

END DO

RETURN

END SUBROUTINE f_computerh

SUBROUTINE f_computeabsvort(u,v,msfu,msfv,msft,cor,dx,dy,av,nx,ny,nz,nxp1,nyp1)

IMPLICIT NONE

INTEGER, INTENT(IN) :: nx,ny,nz,nxp1,nyp1

REAL(KIND=8),DIMENSION(nxp1,ny,nz),INTENT(IN) :: u

REAL(KIND=8),DIMENSION(nx,nyp1,nz),INTENT(IN) :: v

REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(OUT) :: av

REAL(KIND=8),DIMENSION(nxp1,ny),INTENT(IN):: msfu

REAL(KIND=8),DIMENSION(nx,nyp1),INTENT(IN) :: msfv

REAL(KIND=8),DIMENSION(nx,ny),INTENT(IN) :: msft

REAL(KIND=8),DIMENSION(nx,ny),INTENT(IN) :: cor

REAL(KIND=8) :: dx,dy

INTEGER :: jp1,jm1,ip1,im1,i,j,k

REAL(KIND=8) :: dsy,dsx,dudy,dvdx,avort

REAL(KIND=8) :: mm

! PRINT*,'nx,ny,nz,nxp1,nyp1'

! 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)

! PRINT *,jp1,jm1,ip1,im1

dsx = (ip1-im1)*dx

dsy = (jp1-jm1)*dy

mm = msft(i,j)*msft(i,j)

! 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 SUBROUTINE f_computeabsvort

SUBROUTINE f_computepvo(u,v,theta,prs,msfu,msfv,msft,cor,dx,dy,pv,nx,ny,nz,nxp1,nyp1)

IMPLICIT NONE

INTEGER,INTENT(IN) :: nx,ny,nz,nxp1,nyp1

REAL(KIND=8),DIMENSION(nxp1,ny,nz),INTENT(IN) :: u

REAL(KIND=8),DIMENSION(nx,nyp1,nz),INTENT(IN) :: v

REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(IN) :: prs

REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(IN) :: theta

REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(OUT) :: pv

REAL(KIND=8),DIMENSION(nxp1,ny),INTENT(IN) :: msfu

REAL(KIND=8),DIMENSION(nx,nyp1),INTENT(IN) :: msfv

REAL(KIND=8),DIMENSION(nx,ny),INTENT(IN) :: msft

REAL(KIND=8),DIMENSION(nx,ny),INTENT(IN) :: cor

REAL(KIND=8) :: dx,dy

INTEGER :: kp1,km1,jp1,jm1,ip1,im1,i,j,k

REAL(KIND=8) :: dsy,dsx,dp,dudy,dvdx,dudp,dvdp,dthdp,avort

REAL(KIND=8) :: dthdx,dthdy,mm

! PRINT*,'nx,ny,nz,nxp1,nyp1'

! 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)

! PRINT *,jp1,jm1,ip1,im1

dsx = (ip1-im1)*dx

dsy = (jp1-jm1)*dy

mm = msft(i,j)*msft(i,j)

! 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

! if(i.eq.300 .and. j.eq.300) then

! PRINT*,'avort,dudp,dvdp,dthdp,dthdx,dthdy,pv'

! PRINT*,avort,dudp,dvdp,dthdp,dthdx,dthdy,pv(i,j,k)

! endif

pv(i,j,k) = pv(i,j,k)*1.D2

END DO

END DO

END DO

RETURN

END SUBROUTINE f_computepvo

! Theta-e

SUBROUTINE f_computeeth(qvp,tmk,prs,eth,miy,mjx,mkzh)

IMPLICIT NONE

! Input variables

! Qvapor [g/kg]

REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(IN) :: qvp

! Temperature [K]

REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(IN) :: tmk

! full pressure (=P+PB) [hPa]

REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(IN) :: prs

!

! Output variable

! equivalent potential temperature [K]

REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(OUT) :: eth

! Sizes

INTEGER,INTENT(IN) :: miy,mjx,mkzh

! local variables

REAL(KIND=8) :: q

REAL(KIND=8) :: t

REAL(KIND=8) :: p

REAL(KIND=8) :: e

REAL(KIND=8) :: tlcl

INTEGER :: i,j,k

! parameters

REAL(KIND=8),PARAMETER :: EPS = 0.622D0

REAL(KIND=8),PARAMETER :: RGAS = 287.04D0

REAL(KIND=8),PARAMETER :: RGASMD = .608D0

REAL(KIND=8),PARAMETER :: CP = 1004.D0

REAL(KIND=8),PARAMETER :: CPMD = .887D0

REAL(KIND=8),PARAMETER :: GAMMA = RGAS/CP

REAL(KIND=8),PARAMETER :: GAMMAMD = RGASMD - CPMD

REAL(KIND=8),PARAMETER :: TLCLC1 = 2840.D0

REAL(KIND=8),PARAMETER :: TLCLC2 = 3.5D0

REAL(KIND=8),PARAMETER :: TLCLC3 = 4.805D0

REAL(KIND=8),PARAMETER :: TLCLC4 = 55.D0

REAL(KIND=8),PARAMETER :: THTECON1 = 3376.D0

REAL(KIND=8),PARAMETER :: THTECON2 = 2.54D0

REAL(KIND=8),PARAMETER :: THTECON3 = .81D0

DO k = 1,mkzh

DO j = 1,mjx

DO 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))

END DO

END DO

END DO

RETURN

END SUBROUTINE f_computeeth

SUBROUTINE f_computeuvmet(u,v,longca,longcb,flong,flat,&

cen_long,cone,rpd,istag,is_msg_val,umsg,vmsg,uvmetmsg,&

uvmet,nx,ny,nxp1,nyp1,nz)

IMPLICIT NONE

! ISTAG should be 0 if the U,V grids are not staggered.

! That is, NY = NYP1 and NX = NXP1.

INTEGER,INTENT(IN) :: nx,ny,nz,nxp1,nyp1,istag

LOGICAL,INTENT(IN) :: is_msg_val

REAL(KIND=8),DIMENSION(nxp1,ny,nz),INTENT(IN):: u

REAL(KIND=8),DIMENSION(nx,nyp1,nz),INTENT(IN) :: v

REAL(KIND=8),DIMENSION(nx,ny),INTENT(IN) :: flong

REAL(KIND=8),DIMENSION(nx,ny),INTENT(IN) :: flat

REAL(KIND=8),DIMENSION(nx,ny),INTENT(INOUT) :: longca

REAL(KIND=8),DIMENSION(nx,ny),INTENT(INOUT) :: longcb

REAL(KIND=8),INTENT(IN) :: cen_long,cone,rpd

REAL(KIND=8),INTENT(IN) :: umsg,vmsg,uvmetmsg

REAL(KIND=8),DIMENSION(nx,ny,nz,2),INTENT(OUT) :: uvmet

INTEGER :: i,j,k

REAL(KIND=8) :: uk,vk

! msg stands for missing value in this code

! 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

! WRITE (6,FMT=*) ' computing velocities '

DO k = 1,nz

DO j = 1,ny

DO i = 1,nx

IF (istag.EQ.1) THEN

IF (is_msg_val .AND. (u(i,j,k) .EQ. umsg .OR. v(i,j,k) .EQ. vmsg &

.OR. u(i+1,j,k) .EQ. umsg .OR. v(i,j+1,k) .EQ. vmsg)) THEN

uvmet(i,j,k,1) = uvmetmsg

uvmet(i,j,k,2) = uvmetmsg

ELSE

uk = 0.5D0* (u(i,j,k)+u(i+1,j,k))

vk = 0.5D0* (v(i,j,k)+v(i,j+1,k))

uvmet(i,j,k,1) = vk*longcb(i,j) + uk*longca(i,j)

uvmet(i,j,k,2) = vk*longca(i,j) - uk*longcb(i,j)

END IF

ELSE

IF (is_msg_val .AND. (u(i,j,k) .EQ. umsg .OR. v(i,j,k) .EQ. vmsg)) THEN

uvmet(i,j,k,1) = uvmetmsg

uvmet(i,j,k,2) = uvmetmsg

ELSE

uk = u(i,j,k)

vk = v(i,j,k)

uvmet(i,j,k,1) = vk*longcb(i,j) + uk*longca(i,j)

uvmet(i,j,k,2) = vk*longca(i,j) - uk*longcb(i,j)

END IF

END IF

END DO

END DO

END DO

RETURN

END SUBROUTINE f_computeuvmet

SUBROUTINE f_computeomega(qvp,tmk,www,prs,omg,mx,my,mz)

IMPLICIT NONE

INTEGER,INTENT(IN) :: mx, my, mz

REAL(KIND=8),INTENT(IN),DIMENSION(mx,my,mz) :: qvp

REAL(KIND=8),INTENT(IN),DIMENSION(mx,my,mz) :: tmk

REAL(KIND=8),INTENT(IN),DIMENSION(mx,my,mz) :: www

REAL(KIND=8),INTENT(IN),DIMENSION(mx,my,mz) :: prs

REAL(KIND=8),INTENT(OUT),DIMENSION(mx,my,mz) :: omg

! Local variables

INTEGER :: i, j, k

REAL(KIND=8),PARAMETER :: GRAV=9.81, RGAS=287.04, EPS=0.622

DO k=1,mz

DO j=1,my

DO i=1,mx

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)

END DO

END DO

END DO

!

RETURN

END SUBROUTINE f_computeomega

SUBROUTINE f_computetv(temp,qv,tv,nx,ny,nz)

IMPLICIT NONE

INTEGER, INTENT(IN) :: nx,ny,nz

REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(IN) :: temp

REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(IN) :: qv

REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(OUT) :: tv

INTEGER :: i,j,k

REAL(KIND=8),PARAMETER :: EPS = 0.622D0

DO k=1,nz

DO j=1,ny

DO i=1,nx

tv(i,j,k) = temp(i,j,k) * (EPS+qv(i,j,k)) / (EPS * (1.D0+qv(i,j,k)))

END DO

END DO

END DO

RETURN

END SUBROUTINE f_computetv

! Need to deal with the fortran stop statements

SUBROUTINE f_computewetbulb(prs,tmk,qvp,PSADITHTE,PSADIPRS,PSADITMK,throw_exception,twb,nx,ny,nz)

IMPLICIT NONE

EXTERNAL throw_exception

INTEGER,INTENT(IN) :: nx, ny, nz

REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(IN) :: prs

REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(IN) :: tmk

REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(IN) :: qvp

REAL(KIND=8),DIMENSION(150),INTENT(IN) :: PSADITHTE

REAL(KIND=8),DIMENSION(150),INTENT(IN) ::PSADIPRS

REAL(KIND=8),DIMENSION(150,150),INTENT(IN) :: PSADITMK

REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(OUT) :: twb

!EXTERNAL throw_exception

INTEGER :: i,j,k

INTEGER :: jtch,jt,ipch,ip

REAL(KIND=8) :: q, t, p, e, tlcl, eth

REAL(KIND=8) :: fracip,fracip2,fracjt,fracjt2

REAL(KIND=8) :: tonpsadiabat

REAL(KIND=8),PARAMETER :: EPS=0.622

REAL(KIND=8),PARAMETER :: RGAS=287.04

REAL(KIND=8),PARAMETER :: RGASMD=.608

REAL(KIND=8),PARAMETER :: CP=1004.

REAL(KIND=8),PARAMETER :: CPMD=.887

REAL(KIND=8),PARAMETER :: GAMMA=RGAS/CP

REAL(KIND=8),PARAMETER :: GAMMAMD=RGASMD-CPMD

REAL(KIND=8),PARAMETER :: CELKEL=273.15

REAL(KIND=8),PARAMETER :: TLCLC1=2840.

REAL(KIND=8),PARAMETER :: TLCLC2=3.5

REAL(KIND=8),PARAMETER :: TLCLC3=4.805

REAL(KIND=8),PARAMETER :: TLCLC4=55.

REAL(KIND=8),PARAMETER :: THTECON1=3376.

REAL(KIND=8),PARAMETER :: THTECON2=2.54

REAL(KIND=8),PARAMETER :: THTECON3=.81

DO k=1,nz

DO j=1,ny

DO i=1,nx

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))

! Now we need to find the temperature (in K) on a moist adiabat

! (specified by eth in K) given pressure in hPa. It uses a

! lookup table, with data that was generated by the Bolton (1980)

! formula for theta_e.

! First check if pressure is less than min pressure in lookup table.

! If it is, assume parcel is so dry that the given theta-e value can

! be interpretted as theta, and get temperature from the simple dry

! theta formula.

!

IF (p.LE.psadiprs(150)) THEN

tonpsadiabat=eth*(p/1000.)**GAMMA

ELSE

! Otherwise, look for the given thte/prs point in the lookup table.

jt=-1

DO jtch=1,150-1

IF (eth.GE.psadithte(jtch).AND.eth.LT.psadithte(jtch+1)) THEN

jt=jtch

EXIT

ENDIF

END DO

! jt=-1

! 213 CONTINUE

ip=-1

DO ipch=1,150-1

IF (p.LE.psadiprs(ipch).AND.p.GT.psadiprs(ipch+1)) THEN

ip=ipch

EXIT

ENDIF

END DO

! ip=-1

! 215 CONTINUE

IF (jt.EQ.-1.OR.ip.EQ.-1) THEN

CALL throw_exception('Outside of lookup table bounds. prs,thte=',p,eth)

!STOP ! TODO: Need to make python throw an exception here

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.'

CALL throw_exception('Prs and Thte probably unreasonable. prs,thte=',p,eth)

!STOP ! TODO: Need to make python throw an exception here

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

END DO

END DO

END DO

RETURN

END SUBROUTINE f_computewetbulb

SUBROUTINE f_computesrh(u, v, ght, ter, top, sreh, miy, mjx, mkzh)

IMPLICIT NONE

INTEGER,INTENT(IN) :: miy, mjx, mkzh

REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(IN) :: u, v, ght

REAL(KIND=8),INTENT(IN) :: top

REAL(KIND=8),DIMENSION(miy,mjx),INTENT(IN) :: ter

REAL(KIND=8),DIMENSION(miy,mjx),INTENT(OUT) :: sreh

! This helicity code was provided by Dr. Craig Mattocks, and

! verified by Cindy Bruyere to produce results equivalent to

! those generated by RIP4. (The code came from RIP4?)

REAL(KIND=8) :: dh, sdh, su, sv, ua, va, asp, adr, bsp, bdr

REAL(KIND=8) :: cu, cv, x, sum

INTEGER :: i, j, k, k10, k3, ktop

REAL(KIND=8),PARAMETER :: PI=3.14159265d0, DTR=PI/180.d0, DPR=180.d0/PI

DO j = 1, mjx-1

DO i = 1, miy-1

sdh = 0.d0

su = 0.d0

sv = 0.d0

k3 = 0

k10 = 0

ktop = 0

DO k = mkzh, 2, -1

IF (((ght(i,j,k) - ter(i,j)) .GT. 10000.d0) .AND.(k10 .EQ. 0)) THEN

k10 = k

EXIT

ENDIF

IF (((ght(i,j,k) - ter(i,j)) .GT. top) .AND.(ktop .EQ. 0)) THEN

ktop = k

ENDIF

IF (((ght(i,j,k) - ter(i,j)) .GT. 3000.d0) .AND.(k3 .EQ. 0)) THEN

k3 = k

ENDIF

END DO

IF (k10 .EQ. 0) THEN

k10 = 2

ENDIF

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))

END DO

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) THEN

bdr = bdr-360.d0

ENDIF

cu = -bsp * SIN(bdr*dtr)

cv = -bsp * COS(bdr*dtr)

sum = 0.d0

DO 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

END DO

sreh(i,j) = -sum

END DO

END DO

RETURN

END SUBROUTINE f_computesrh

SUBROUTINE f_computeuh(zp,mapfct,dx,dy,uhmnhgt,uhmxhgt,us,vs,w,tem1,tem2,uh,nx,ny,nz,nzp1)

IMPLICIT NONE

INTEGER, INTENT(IN) :: nx,ny,nz,nzp1

REAL(KIND=8),DIMENSION(nx,ny,nzp1),INTENT(IN) :: zp

REAL(KIND=8),DIMENSION(nx,ny),INTENT(IN) :: mapfct

REAL(KIND=8),INTENT(IN) :: dx,dy

REAL(KIND=8),INTENT(IN) :: uhmnhgt,uhmxhgt

REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(IN) :: us

REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(IN) :: vs

REAL(KIND=8),DIMENSION(nx,ny,nzp1),INTENT(IN) :: w

REAL(KIND=8),DIMENSION(nx,ny),INTENT(OUT) :: uh

REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(INOUT) :: tem1

REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(INOUT) :: tem2

!

! Misc local variables

!

INTEGER :: i,j,k,kbot,ktop

REAL(KIND=8) :: twodx,twody,wgtlw,sum,wmean,wsum,wavg

REAL(KIND=8) :: helbot,heltop,wbot,wtop

REAL(KIND=8) :: 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)