/
wrf_fctt.f90
131 lines (110 loc) · 4.21 KB
/
wrf_fctt.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
!NCLFORTSTART
SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, pf, haveqci,&
fill_nocloud, missing, opt_thresh, nz, ns, ew)
USE wrf_constants, ONLY : EPS, USSALR, RD, G, ABSCOEFI, ABSCOEF, CELKEL
IMPLICIT NONE
!f2py threadsafe
!f2py intent(in,out) :: ctt
INTEGER, INTENT(IN) :: nz, ns, ew, haveqci, fill_nocloud
REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(IN) :: ght, prs, tk, qci, qcw, qvp
REAL(KIND=8), DIMENSION(ew,ns), INTENT(IN) :: ter
REAL(KIND=8), DIMENSION(ew,ns), INTENT(OUT) :: ctt
REAL(KIND=8), INTENT(IN) :: missing
REAL(KIND=8), INTENT(IN) :: opt_thresh
REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(INOUT) :: pf
!NCLEND
! REAL(KIND=8) :: znfac(nz)
! LOCAL VARIABLES
INTEGER i,j,k,ripk
REAL(KIND=8) :: opdepthu, opdepthd, dp, arg1, fac, prsctt, ratmix
REAL(KIND=8) :: arg2, agl_hgt, vt
REAL(KIND=8) :: p1, p2
!$OMP PARALLEL
! Calculate the surface pressure
!$OMP DO COLLAPSE(2) PRIVATE(i, j, ratmix, arg1, arg2, vt, agl_hgt) &
!$OMP SCHEDULE(runtime)
DO j=1,ns
DO i=1,ew
ratmix = .001D0*qvp(i,j,1)
arg1 = EPS + ratmix
arg2 = EPS*(1. + ratmix)
vt = tk(i,j,1)*arg1/arg2 !Virtual temperature
agl_hgt = ght(i,j,nz) - ter(i,j)
arg1 = -G/(RD*USSALR)
pf(i,j,nz) = prs(i,j,1)*(vt/(vt + USSALR*(agl_hgt)))**(arg1)
END DO
END DO
!$OMP END DO
!$OMP DO COLLAPSE(3) PRIVATE(i, j, k, ripk) SCHEDULE(runtime)
DO k=1,nz-1
DO j=1,ns
DO i=1,ew
ripk = nz-k+1
pf(i,j,k) = .5D0*(prs(i,j,ripk) + prs(i,j,ripk-1))
END DO
END DO
END DO
!$OMP END DO
!$OMP DO COLLAPSE(2) PRIVATE(i, j, k, ripk, opdepthd, opdepthu, &
!$OMP prsctt, dp, p1, p2, fac, arg1) SCHEDULE(runtime)
DO j=1,ns
DO i=1,ew
opdepthd = 0.D0
k = 0
prsctt = -1
! Integrate downward from model top, calculating path at full
! model vertical levels.
DO k=2,nz
opdepthu = opdepthd
ripk = nz - k + 1
IF (k .NE. 1) THEN
dp = 100.D0*(pf(i,j,k) - pf(i,j,k-1)) ! should be in Pa
ELSE
dp = 200.D0*(pf(i,j,1) - prs(i,j,nz)) ! should be in Pa
END IF
IF (haveqci .EQ. 0) then
IF (tk(i,j,ripk) .LT. CELKEL) then
! Note: abscoefi is m**2/g, qcw is g/kg, so no convrsion needed
opdepthd = opdepthu + ABSCOEFI*qcw(i,j,ripk) * dp/G
ELSE
opdepthd = opdepthu + ABSCOEF*qcw(i,j,ripk) * dp/G
END IF
ELSE
opdepthd = opdepthd + (ABSCOEF*qcw(i,j,ripk) + ABSCOEFI*qci(i,j,ripk))*dp/G
END IF
IF (opdepthd .LT. opt_thresh .AND. k .LT. nz) THEN
CYCLE
ELSE IF (opdepthd .LT. opt_thresh .AND. k .EQ. nz) THEN
IF (fill_nocloud .EQ. 0) THEN
prsctt = prs(i,j,1)
ENDIF
EXIT
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(i,j,1), MAX(prs(i,j,nz), prsctt))
EXIT
END IF
END DO
! prsctt should only be 0 if fill values are used
IF (prsctt .GT. -1) THEN
DO k=2,nz
ripk = nz - k + 1
p1 = prs(i,j,ripk+1)
p2 = prs(i,j,ripk)
IF (prsctt .GE. p1 .AND. prsctt .LE. p2) THEN
fac = (prsctt - p1)/(p2 - p1)
arg1 = fac*(tk(i,j,ripk) - tk(i,j,ripk+1)) - CELKEL
ctt(i,j) = tk(i,j,ripk+1) + arg1
EXIT
END IF
END DO
ELSE
ctt(i,j) = missing
END IF
END DO
END DO
!$OMP END DO
!$OMP END PARALLEL
RETURN
END SUBROUTINE wrfcttcalc