-
Notifications
You must be signed in to change notification settings - Fork 237
/
fizhi_rayleigh.F
130 lines (115 loc) · 3.5 KB
/
fizhi_rayleigh.F
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
#include "FIZHI_OPTIONS.h"
subroutine rayleigh(myid,pres,pkap,pekap,zsurf,u,v,t,s,im,jm,lm,
. bi,bj,rfu,rfv,rft)
C **********************************************************************
C
C PURPOSE
C Rayleigh Friction -- Linear Drag (Strong Damping) Above 70 Km
C
C ARGUMENTS DESCRIPTION
C
C INPUT:
C MYID .... PROCESS(OR) NUMBER
C PRES .... MID-LEVEL PRESSURE IN MB
C PKAP .... MID-LEVEL PRESSURE ** KAPPA
C PEKAP ... EDGE-LEVEL PRESSURE ** KAPPA
C ZSURF ... SURFACE ELEVATION IN M
C U ....... U-WIND
C V ....... V-WIND
C TH ...... THETA (ACTUALLY REAL THETA * P0**KAPPA) IN K
C S ...... SPECIFIC HUMIDITY (KG/KG)
C IM ...... NUMBER OF LONGITUDE POINTS
C JM ...... NUMBER OF LATITUDE POINTS
C LM ...... NUMBER OF VERTICAL LEVELS
C BI ...... X-DIRECTION PROCESSOR INDEX
C BJ ...... Y-DIRECTION PROCESSOR INDEX
C OUTPUT:
C RFU ..... U-WIND TENDENCY
C RFV ..... V-WIND TENDENCY
C RFT ..... THETA TENDENCY
C
C **********************************************************************
implicit none
integer myid,im,jm,lm,bi,bj
_RL zsurf(im,jm),pres(im,jm,lm),pkap(im,jm,lm)
_RL pekap(im,jm,lm+1)
_RL u(im,jm,lm),v(im,jm,lm),t(im,jm,lm),s(im,jm,lm)
_RL rfu(im,jm,lm),rfv(im,jm,lm),rft(im,jm,lm)
integer i,j,L
_RL rf(im,jm,lm)
_RL z(im,jm,lm)
_RL dz(im,jm,lm)
_RL cpog, cpinv, virtcon, getcon, dampcoef
#ifdef ALLOW_DIAGNOSTICS
logical diagnostics_is_on
external diagnostics_is_on
_RL tmpdiag(im,jm)
#endif
C **********************************************************************
C **** APPLY RAYLEIGH FRICTION TO WIND (INCLUDE HEATING) ***
C **********************************************************************
cpog = getcon('CP')/getcon('GRAVITY')
cpinv = 1.0/getcon('CP')
virtcon = getcon('VIRTCON')
dampcoef = 2./3.
do L=1,lm
do j=1,jm
do i=1,im
dz(i,j,L) = cpog * (pekap(i,j,L+1)-pekap(i,j,L)) * t(i,j,L) *
. (1.+virtcon*s(i,j,L))
enddo
enddo
enddo
do j=1,jm
do i=1,im
z(i,j,lm) = zsurf(i,j) + 0.5 * dz(i,j,lm)
enddo
enddo
do L=lm-1,1,-1
do j=1,jm
do i=1,im
z(i,j,L) = z(i,j,L+1) + 0.5 * (dz(i,j,L)+dz(i,j,L+1))
enddo
enddo
enddo
do L=1,lm
do j=1,jm
do i=1,im
rf(i,j,L) = dampcoef*(1+tanh((z(i,j,L)-50000.)/5000.))/86400.
rfu(i,j,L) = - rf(i,j,L) * u(i,j,L)
rfv(i,j,L) = - rf(i,j,L) * v(i,j,L)
rft(i,j,L) = -(u(i,j,L)*rfu(i,j,L) + v(i,j,L)*rfv(i,j,L) )*cpinv
. /pkap(i,j,L)
enddo
enddo
enddo
#ifdef ALLOW_DIAGNOSTICS
do L=1,lm
if(diagnostics_is_on('RFU ',myid) ) then
do j=1,jm
do i=1,im
tmpdiag(i,j) = rfu(i,j,L)*86400
enddo
enddo
C call diagnostics_fill(tmpdiag,'RFU ',L,1,3,bi,bj,myid)
endif
if(diagnostics_is_on('RFV ',myid) ) then
do j=1,jm
do i=1,im
tmpdiag(i,j) = rfv(i,j,L)*86400
enddo
enddo
C call diagnostics_fill(tmpdiag,'RFV ',L,1,3,bi,bj,myid)
endif
if(diagnostics_is_on('RFT ',myid) ) then
do j=1,jm
do i=1,im
tmpdiag(i,j) = rft(i,j,L)*86400
enddo
enddo
C call diagnostics_fill(tmpdiag,'RFT ',L,1,3,bi,bj,myid)
endif
enddo
#endif
return
end