forked from 3rdparty/wrf-python
You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
956 lines
26 KiB
956 lines
26 KiB
! NCLFORTSTART |
|
SUBROUTINE DCOMPUTEPI(pi, pressure, nx, ny, nz) |
|
USE wrf_constants, ONLY : P1000MB, RD, CP |
|
|
|
IMPLICIT NONE |
|
|
|
!f2py threadsafe |
|
!f2py intent(in,out) :: pi |
|
|
|
INTEGER, INTENT(IN) :: nx, ny, nz |
|
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(OUT) :: pi |
|
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: pressure |
|
! NCLEND |
|
|
|
INTEGER i, j, k |
|
|
|
!REAL(KIND=8), PARAMETER :: P1000MB=100000.D0, R_D=287.D0, CP=7.D0*R_D/2.D0 |
|
|
|
!$OMP PARALLEL DO COLLAPSE(3) |
|
DO k = 1,nz |
|
DO j = 1,ny |
|
DO i = 1,nx |
|
pi(i,j,k) = (pressure(i,j,k)/P1000MB)**(RD/CP) |
|
END DO |
|
END DO |
|
END DO |
|
!$OMP END PARALLEL DO |
|
|
|
END SUBROUTINE DCOMPUTEPI |
|
|
|
! Temperature from potential temperature in kelvin. |
|
!NCLFORTSTART |
|
SUBROUTINE DCOMPUTETK(tk, pressure, theta, nx) |
|
USE wrf_constants, ONLY : P1000MB, RD, CP |
|
|
|
IMPLICIT NONE |
|
|
|
!f2py threadsafe |
|
!f2py intent(in,out) :: tk |
|
|
|
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 |
|
|
|
! NCLEND |
|
|
|
INTEGER :: i |
|
|
|
!REAL(KIND=8), PARAMETER :: P1000MB=100000.D0, RD=287.D0, CP=7.D0*RD/2.D0 |
|
|
|
!$OMP PARALLEL DO |
|
DO i = 1,nx |
|
!pi = (pressure(i)/P1000MB)**(RD/CP) |
|
!tk(i) = pi * theta(i) |
|
tk(i) = (pressure(i)/P1000MB)**(RD/CP) * theta(i) |
|
END DO |
|
!$OMP END PARALLEL DO |
|
|
|
RETURN |
|
|
|
END SUBROUTINE DCOMPUTETK |
|
|
|
|
|
! NCLFORTSTART |
|
SUBROUTINE DINTERP3DZ(data3d, out2d, zdata, desiredloc, nx, ny, nz, missingval) |
|
IMPLICIT NONE |
|
|
|
!f2py threadsafe |
|
!f2py intent(in,out) :: out2d |
|
|
|
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 |
|
|
|
! NCLEND |
|
|
|
INTEGER :: i,j,kp,ip,im |
|
LOGICAL :: dointerp |
|
REAL(KIND=8) :: w1,w2 |
|
|
|
! 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 |
|
|
|
!$OMP PARALLEL DO COLLAPSE(2) PRIVATE(i,j,kp,dointerp,w1,w2) & |
|
!$OMP FIRSTPRIVATE(ip,im) |
|
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) < desiredloc) .AND. (zdata(i,j,kp-ip) > desiredloc))) THEN |
|
w2 = (desiredloc - 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 |
|
!$OMP END PARALLEL DO |
|
|
|
RETURN |
|
|
|
END SUBROUTINE DINTERP3DZ |
|
|
|
! PORT DZSTAG HERE |
|
|
|
! NCLFORTSTART |
|
SUBROUTINE DZSTAG(znew, nx, ny, nz, z, nxz, nyz ,nzz, terrain) |
|
IMPLICIT NONE |
|
|
|
!f2py threadsafe |
|
!f2py intent(in,out) :: znew |
|
|
|
INTEGER, INTENT(IN) :: nx, ny, nz, nxz, nyz, nzz |
|
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(OUT) :: znew |
|
REAL(KIND=8), DIMENSION(nxz,nyz,nzz), INTENT(IN) :: z |
|
REAL(KIND=8), DIMENSION(nxz,nyz), INTENT(IN) :: terrain |
|
! NCLEND |
|
|
|
INTEGER :: i,j,k,ii,im1,jj,jm1 |
|
|
|
! check for u, v, or w (x,y,or z) staggering |
|
! for x and y stag, avg z to x, y, point |
|
IF (nx .GT. nxz) THEN |
|
DO k = 1,nz |
|
DO j = 1,ny |
|
DO i = 1,nx |
|
ii = MIN(i,nxz) |
|
im1 = MAX(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 = MIN(j,nyz) |
|
jm1 = MAX(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 |
|
|
|
! w (z) staggering |
|
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 SUBROUTINE DZSTAG |
|
|
|
! NCLFORTSTART |
|
SUBROUTINE DINTERP2DXY(v3d, v2d, xy, nx, ny, nz, nxy) |
|
|
|
IMPLICIT NONE |
|
|
|
!f2py threadsafe |
|
!f2py intent(in,out) :: v2d |
|
|
|
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 |
|
|
|
! NCLEND |
|
|
|
INTEGER :: i, j, k, ij |
|
REAL(KIND=8) :: w11, w12, w21, w22, wx, wy |
|
|
|
!$OMP PARALLEL DO PRIVATE(i,j,k,ij,w11,w12,w21,w22,wx,wy) |
|
DO ij = 1,nxy |
|
i = MAX(1,MIN(nx-1,INT(xy(1,ij)+1))) |
|
j = MAX(1,MIN(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 |
|
!$OMP END PARALLEL DO |
|
|
|
RETURN |
|
|
|
END SUBROUTINE DINTERP2DXY |
|
|
|
|
|
! NCLFORTSTART |
|
SUBROUTINE DINTERP1D(v_in, v_out, z_in, z_out, vmsg, nz_in, nz_out) |
|
|
|
IMPLICIT NONE |
|
|
|
!f2py threadsafe |
|
!f2py intent(in,out) :: v_out |
|
|
|
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 |
|
|
|
! NCLEND |
|
|
|
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 |
|
|
|
!$OMP PARALLEL DO PRIVATE(kp, k, interp, height, w1, w2) FIRSTPRIVATE(ip, im) |
|
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 |
|
!$OMP END PARALLEL DO |
|
|
|
RETURN |
|
|
|
END SUBROUTINE DINTERP1D |
|
|
|
! 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 |
|
! |
|
|
|
! NCLFORTSTART |
|
SUBROUTINE DCOMPUTESEAPRS(nx, ny, nz, z, t, p, q, sea_level_pressure, & |
|
t_sea_level, t_surf, level, errstat, errmsg) |
|
USE wrf_constants, ONLY : ALGERR, RD, G, USSALR |
|
|
|
IMPLICIT NONE |
|
|
|
!f2py threadsafe |
|
!f2py intent(in,out) :: sea_level_pressure |
|
|
|
! 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 |
|
INTEGER, INTENT(INOUT) :: errstat |
|
CHARACTER(LEN=*), INTENT(INOUT) :: errmsg |
|
|
|
! NCLFORTEND |
|
|
|
! 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 |
|
INTEGER :: errcnt, bad_i, bad_j, bad_sfp |
|
|
|
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. |
|
|
|
errstat = 0 |
|
errcnt = 0 |
|
bad_i = -1 |
|
bad_j = -1 |
|
bad_sfp = -1 |
|
|
|
!$OMP PARALLEL DO COLLAPSE(2) PRIVATE(i,j,k,found) REDUCTION(+:errcnt) |
|
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 |
|
errcnt = errcnt + 1 |
|
!$OMP CRITICAL |
|
! Only do this the first time |
|
IF (bad_i .EQ. -1) THEN |
|
bad_i = i |
|
bad_j = j |
|
bad_sfp = p(i,j,1) / 100. |
|
END IF |
|
!$OMP END CRITICAL |
|
END IF |
|
END DO |
|
END DO |
|
!$OMP END PARALLEL DO |
|
|
|
IF (errcnt > 0) THEN |
|
errstat = ALGERR |
|
WRITE(errmsg,*) "Error in finding 100 hPa up. i=", bad_i, "j=", bad_j, "sfc_p=", bad_sfp |
|
RETURN |
|
END IF |
|
|
|
! Get temperature PCONST Pa above surface. Use this to extrapolate |
|
! the temperature at the surface and down to sea level. |
|
!$OMP PARALLEL DO COLLAPSE(2) PRIVATE(i,j,klo,khi) REDUCTION(+:errcnt) |
|
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 |
|
errcnt = errcnt + 1 |
|
!$OMP CRITICAL |
|
IF (bad_i .EQ. -1) THEN |
|
bad_i = i |
|
bad_j = j |
|
END IF |
|
!$OMP END CRITICAL |
|
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 = 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)**(USSALR*RD/G) |
|
t_sea_level(i,j) = t_at_pconst + USSALR*z_at_pconst |
|
END DO |
|
END DO |
|
!$OMP END PARALLEL DO |
|
|
|
IF (errcnt > 0) THEN |
|
errstat = ALGERR |
|
WRITE(errmsg,*) "Error trapping levels at i=", bad_i, "j=", bad_j |
|
RETURN |
|
END IF |
|
|
|
! 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 |
|
!$OMP PARALLEL DO COLLAPSE(2) PRIVATE(l1,l2,l3) |
|
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 |
|
!$OMP END PARALLEL DO |
|
END IF |
|
|
|
! The grand finale: ta da! |
|
!$OMP PARALLEL DO COLLAPSE(2) |
|
DO j = 1,ny |
|
DO i = 1,nx |
|
!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(i,j,1))/& |
|
(RD*(t_sea_level(i,j) + t_surf(i,j))))) |
|
END DO |
|
END DO |
|
!$OMP END PARALLEL 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 DCOMPUTESEAPRS |
|
|
|
|
|
! Double precision version. If you make a change here, you |
|
! must make the same change below to filter2d. |
|
|
|
! NCLFORTSTART |
|
SUBROUTINE DFILTER2D(a, b, nx, ny, it, missing) |
|
|
|
IMPLICIT NONE |
|
|
|
!f2py threadsafe |
|
!f2py intent(in,out) :: a |
|
|
|
INTEGER, INTENT(IN) :: nx, ny, it |
|
REAL(KIND=8), DIMENSION(nx, ny), INTENT(INOUT) :: a |
|
REAL(KIND=8), INTENT(IN) :: missing |
|
REAL(KIND=8), DIMENSION(nx, ny), INTENT(INOUT) :: b |
|
|
|
! NCLEND |
|
|
|
REAL(KIND=8), PARAMETER :: COEF=0.25D0 |
|
|
|
INTEGER :: i, j, iter |
|
|
|
DO iter=1,it |
|
!$OMP PARALLEL DO COLLAPSE(2) |
|
DO j=1,ny |
|
DO i = 1,nx |
|
b(i,j) = a(i,j) |
|
END DO |
|
END DO |
|
!$OMP END PARALLEL DO |
|
|
|
!$OMP PARALLEL DO COLLAPSE(2) |
|
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 |
|
!$OMP END PARALLEL DO |
|
|
|
!$OMP PARALLEL DO COLLAPSE(2) |
|
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 |
|
!$OMP END PARALLEL DO |
|
! 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 |
|
! a(i,j) = a(i,j) - .99*coef*(b(i,j-1)-2*b(i,j)+b(i,j+1)) |
|
! enddo |
|
! enddo |
|
! do j=1,ny |
|
! do i=2,nx-1 |
|
! a(i,j) = a(i,j) - .99*coef*(b(i-1,j)-2*b(i,j)+b(i+1,j)) |
|
! enddo |
|
! enddo |
|
END DO |
|
|
|
RETURN |
|
|
|
END SUBROUTINE DFILTER2D |
|
|
|
! Single precision version. If you make a change here, you |
|
! must make the same change below to dfilter2d. |
|
|
|
! NCLFORTSTART |
|
SUBROUTINE FILTER2D(a, b, nx, ny, it, missing) |
|
IMPLICIT NONE |
|
|
|
!f2py threadsafe |
|
!f2py intent(in,out) :: a |
|
|
|
INTEGER, INTENT(IN) :: nx, ny, it |
|
REAL(KIND=4), DIMENSION(nx, ny), INTENT(INOUT) :: a |
|
REAL(KIND=4), INTENT(IN) :: missing |
|
REAL(KIND=4), DIMENSION(nx, ny), INTENT(INOUT) :: b |
|
|
|
! NCLEND |
|
|
|
REAL(KIND=4), PARAMETER :: COEF=0.25 |
|
|
|
INTEGER :: i, j, iter |
|
|
|
!$OMP PARALLEL |
|
|
|
DO iter=1,it |
|
!$OMP DO COLLAPSE(2) |
|
DO j=1,ny |
|
DO i = 1,nx |
|
b(i,j) = a(i,j) |
|
END DO |
|
END DO |
|
!$OMP END DO |
|
|
|
!$OMP DO COLLAPSE(2) |
|
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 |
|
!$OMP END DO |
|
|
|
!$OMP DO COLLAPSE(2) |
|
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 |
|
!$OMP END DO |
|
! 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 |
|
! a(i,j) = a(i,j) - .99*coef*(b(i,j-1)-2*b(i,j)+b(i,j+1)) |
|
! enddo |
|
! enddo |
|
! do j=1,ny |
|
! do i=2,nx-1 |
|
! a(i,j) = a(i,j) - .99*coef*(b(i-1,j)-2*b(i,j)+b(i+1,j)) |
|
! enddo |
|
! enddo |
|
END DO |
|
|
|
!$OMP END PARALLEL |
|
|
|
RETURN |
|
|
|
END SUBROUTINE FILTER2D |
|
|
|
|
|
! NCLFORTSTART |
|
SUBROUTINE DCOMPUTERH(qv, p, t, rh, nx) |
|
USE wrf_constants, ONLY : EZERO, ESLCON1, ESLCON2, CELKEL, RD, RV, EPS |
|
|
|
IMPLICIT NONE |
|
|
|
!f2py threadsafe |
|
!f2py intent(in,out) :: rh |
|
|
|
INTEGER, INTENT(IN) :: nx |
|
REAL(KIND=8), DIMENSION(nx), INTENT(IN) :: qv, p, t |
|
REAL(KIND=8), DIMENSION(nx), INTENT(OUT) :: rh |
|
|
|
! NCLEND |
|
|
|
INTEGER :: i |
|
REAL(KIND=8) :: qvs,es,pressure,temperature |
|
|
|
!$OMP PARALLEL DO PRIVATE(qvs, es, pressure, temperature) |
|
DO i = 1,nx |
|
pressure = p(i) |
|
temperature = t(i) |
|
! es = 1000.*svp1* |
|
es = EZERO*EXP(ESLCON1*(temperature - CELKEL)/(temperature - ESLCON2)) |
|
! qvs = ep_2*es/(pressure-es) |
|
qvs = EPS*es/(0.01D0*pressure - (1.D0 - EPS)*es) |
|
! rh = 100*MAX(1., qv(i)/qvs) |
|
! rh(i) = 100.*qv(i)/qvs |
|
rh(i) = 100.D0*MAX(MIN(qv(i)/qvs, 1.0D0), 0.0D0) |
|
END DO |
|
!$OMP END PARALLEL DO |
|
|
|
RETURN |
|
|
|
END SUBROUTINE DCOMPUTERH |
|
|
|
|
|
! NCLFORTSTART |
|
SUBROUTINE DGETIJLATLONG(lat_array, long_array, lat, longitude, ii, jj, nx, ny, imsg) |
|
IMPLICIT NONE |
|
|
|
!f2py threadsafe |
|
!f2py intent(in,out) :: ii,jj |
|
|
|
INTEGER, INTENT(IN) :: nx, ny, imsg |
|
INTEGER, INTENT(OUT) :: ii, jj |
|
REAL(KIND=8), DIMENSION(nx,ny), INTENT(IN) :: lat_array, long_array |
|
REAL(KIND=8) :: lat, longitude |
|
! NCLEND |
|
|
|
REAL(KIND=8) :: longd,latd |
|
INTEGER :: i,j |
|
REAL(KIND=8) :: ir,jr |
|
REAL(KIND=8) :: dist_min,dist |
|
|
|
! 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 |
|
! longd = MIN((long_array(i,j)-longitude)**2, & |
|
! (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 |
|
|
|
! The original version of this routine returned IR and JR. But, then |
|
! the NCL script that called this routine was converting IR and JR |
|
! to integer, so why not just return II and JJ? |
|
|
|
! Also, I'm subtracing 1 here, because it will be returned to NCL |
|
! script which has 0-based indexing. |
|
|
|
IF (ir .NE. imsg .AND. jr .NE. imsg) THEN |
|
ii = NINT(ir) - 1 |
|
jj = NINT(jr) - 1 |
|
ELSE |
|
ii = imsg |
|
jj = imsg |
|
END IF |
|
|
|
! we will just return the nearest point at present |
|
|
|
RETURN |
|
|
|
END SUBROUTINE DGETIJLATLONG |
|
|
|
! You need to modify the C-WRAPPER in NCL for this to work |
|
|
|
! 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 |
|
|
|
! ISTAG should be 0 if the U,V grids are not staggered. |
|
! That is, NY = NYP1 and NX = NXP1. |
|
|
|
!f2py threadsafe |
|
!f2py intent(in,out) :: uvmet |
|
|
|
INTEGER,INTENT(IN) :: nx, ny, nxp1, nyp1, istag |
|
LOGICAL,INTENT(IN) :: is_msg_val |
|
REAL(KIND=8), DIMENSION(nxp1,ny), INTENT(IN):: u |
|
REAL(KIND=8), DIMENSION(nx,nyp1), 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,2), INTENT(OUT) :: uvmet |
|
|
|
|
|
! NCLEND |
|
|
|
INTEGER :: i,j |
|
!REAL(KIND=8) :: uk, vk |
|
|
|
! msg stands for missing value in this code |
|
! WRITE (6,FMT=*) ' in compute_uvmet ',NX,NY,NXP1,NYP1,ISTAG |
|
|
|
!$OMP PARALLEL |
|
|
|
!$OMP DO COLLAPSE(2) |
|
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 |
|
!$OMP END DO |
|
|
|
! Note: Intentionally removed as many IF statements as possible from loops |
|
! to improve vectorization. |
|
|
|
IF (istag .EQ. 0) THEN ! Not staggered |
|
IF (.NOT. is_msg_val) THEN ! No missing values used |
|
!$OMP DO COLLAPSE(2) |
|
DO j = 1,ny |
|
DO i = 1,nx |
|
uvmet(i,j,1) = v(i,j)*longcb(i,j) + u(i,j)*longca(i,j) |
|
uvmet(i,j,2) = v(i,j)*longca(i,j) - u(i,j)*longcb(i,j) |
|
END DO |
|
END DO |
|
!$OMP END DO |
|
ELSE ! Missing values used |
|
!$OMP DO COLLAPSE(2) |
|
DO j = 1,ny |
|
DO i = 1,nx |
|
IF ((u(i,j) .NE. umsg .AND. v(i,j) .NE. vmsg)) THEN |
|
uvmet(i,j,1) = v(i,j)*longcb(i,j) + u(i,j)*longca(i,j) |
|
uvmet(i,j,2) = v(i,j)*longca(i,j) - u(i,j)*longcb(i,j) |
|
ELSE |
|
uvmet(i,j,1) = uvmetmsg |
|
uvmet(i,j,2) = uvmetmsg |
|
END IF |
|
END DO |
|
END DO |
|
!$OMP END DO |
|
END IF |
|
ELSE ! Staggered |
|
IF (.NOT. is_msg_val) THEN ! No missing values used |
|
!$OMP DO COLLAPSE(2) |
|
DO j = 1,ny |
|
DO i = 1,nx |
|
! This is the more readable version. |
|
!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) |
|
|
|
uvmet(i,j,1) = (0.5D0*(v(i,j) + v(i,j+1)))*longcb(i,j) + & |
|
(0.5D0*(u(i,j) + u(i+1,j)))*longca(i,j) |
|
uvmet(i,j,2) = (0.5D0*(v(i,j) + v(i,j+1)))*longca(i,j) - & |
|
(0.5D0*(u(i,j) + u(i+1,j)))*longcb(i,j) |
|
|
|
END DO |
|
END DO |
|
!$OMP END DO |
|
ELSE ! Missing values used |
|
!$OMP DO COLLAPSE(2) |
|
DO j = 1,ny |
|
DO i = 1,nx |
|
IF (u(i,j) .NE. umsg .AND. v(i,j) .NE. vmsg .AND. u(i+1,j) .NE. umsg .AND. v(i,j+1) .NE. vmsg) THEN |
|
! This is the more readable version. |
|
!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) |
|
|
|
uvmet(i,j,1) = (0.5D0*(v(i,j) + v(i,j+1)))*longcb(i,j) + & |
|
(0.5D0*(u(i,j) + u(i+1,j)))*longca(i,j) |
|
uvmet(i,j,2) = (0.5D0*(v(i,j) + v(i,j+1)))*longca(i,j) - & |
|
(0.5D0*(u(i,j) + u(i+1,j)))*longcb(i,j) |
|
ELSE |
|
uvmet(i,j,1) = uvmetmsg |
|
uvmet(i,j,2) = uvmetmsg |
|
END IF |
|
END DO |
|
END DO |
|
!$OMP END DO |
|
END IF |
|
END IF |
|
|
|
!$OMP END PARALLEL |
|
|
|
RETURN |
|
|
|
END SUBROUTINE DCOMPUTEUVMET |
|
|
|
|
|
|
|
|
|
! This was originally a routine that took 2D input arrays. Since |
|
! the NCL C wrapper routine can handle multiple dimensions, it's |
|
! not necessary to have anything bigger than 1D here. |
|
|
|
! NCLFORTSTART |
|
SUBROUTINE DCOMPUTETD(td, pressure, qv_in, nx) |
|
IMPLICIT NONE |
|
|
|
!f2py threadsafe |
|
!f2py intent(in,out) :: td |
|
|
|
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 |
|
|
|
! NCLEND |
|
|
|
REAL(KIND=8) :: qv,tdc |
|
|
|
INTEGER :: i |
|
|
|
!$OMP PARALLEL DO PRIVATE(i,qv,tdc) |
|
DO i = 1,nx |
|
qv = MAX(qv_in(i), 0.D0) |
|
! vapor pressure |
|
tdc = qv*pressure(i)/(.622D0 + qv) |
|
|
|
! avoid problems near zero |
|
tdc = MAX(tdc, 0.001D0) |
|
td(i) = (243.5D0*LOG(tdc) - 440.8D0)/(19.48D0 - LOG(tdc)) |
|
END DO |
|
!$OMP END PARALLEL DO |
|
|
|
RETURN |
|
|
|
END SUBROUTINE DCOMPUTETD |
|
|
|
! NCLFORTSTART |
|
SUBROUTINE DCOMPUTEICLW(iclw, pressure, qc_in, nx, ny, nz) |
|
USE wrf_constants, ONLY : G |
|
|
|
IMPLICIT NONE |
|
|
|
!f2py threadsafe |
|
!f2py intent(in,out) :: iclw |
|
|
|
INTEGER, INTENT(IN) :: nx,ny,nz |
|
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: pressure |
|
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: qc_in |
|
REAL(KIND=8), DIMENSION(nx,ny), INTENT(OUT) :: iclw |
|
|
|
! NCLEND |
|
|
|
REAL(KIND=8) :: qclw, dp |
|
REAL(KIND=8), PARAMETER :: GG = 1000.D0/G |
|
INTEGER i,j,k |
|
|
|
iclw = 0 |
|
|
|
DO j = 3,ny - 2 |
|
DO i = 3,nx - 2 |
|
DO k = 1,nz |
|
qclw = MAX(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 SUBROUTINE DCOMPUTEICLW |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|