forked from altMITgcm/MITgcm
-
Notifications
You must be signed in to change notification settings - Fork 0
/
phy_shtorh.F
159 lines (137 loc) · 4.13 KB
/
phy_shtorh.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
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
C $Header: /u/gcmpack/MITgcm/pkg/aim_v23/phy_shtorh.F,v 1.4 2006/06/22 01:25:02 jmc Exp $
C $Name: $
#include "AIM_OPTIONS.h"
SUBROUTINE SHTORH (IMODE,NGP,TA,PS,SIG,QA,RH,QSAT,myThid)
C--
C-- SUBROUTINE SHTORH (IMODE,NGP,TA,PS,SIG,QA,RH,QSAT)
C--
C-- Purpose: compute saturation specific humidity and
C-- relative hum. from specific hum. (or viceversa)
C-- Input: IMODE : mode of operation
C-- NGP : no. of grid-points
C-- TA : abs. temperature
C-- PS : normalized pressure (= p/1000_hPa) [if SIG < 0]
C-- : normalized sfc. pres. (= ps/1000_hPa) [if SIG > 0]
C-- SIG : sigma level
C-- QA : specific humidity in g/kg [if IMODE = 1]
C-- RH : relative humidity [if IMODE < 0]
C-- Output: RH : relative humidity [if IMODE = 1]
C-- QA : specific humidity in g/kg [if IMODE < 0]
C-- QSAT : saturation spec. hum. in g/kg
C-- RH : d.Qsat/d.T in g/kg/K [if IMODE = 2]
C--
IMPLICIT NONE
C-- Routine arguments:
INTEGER IMODE, NGP
INTEGER myThid
c _RL TA(NGP), PS(NGP), QA(NGP), RH(NGP), QSAT(NGP)
_RL TA(NGP), PS(NGP), QSAT(NGP), QA(*), RH(*)
C- jmc: declare all routine arguments:
_RL SIG
#ifdef ALLOW_AIM
C-- Local variables:
INTEGER J
C- jmc: declare all local variables:
_RL E0, C1, C2, T0, T1, T2, QS1, QS2
_RL sigP, recT, tmpQ
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C--- 1. Compute Qsat (g/kg) from T (degK) and normalized pres. P (= p/1000_hPa)
C If SIG > 0, P = Ps * sigma, otherwise P = Ps(1) = const.
C
E0= 6.108 _d -3
C1= 17.269 _d 0
C2= 21.875 _d 0
T0=273.16 _d 0
T1= 35.86 _d 0
T2= 7.66 _d 0
QS1= 622. _d 0
QS2= .378 _d 0
IF (IMODE.EQ.2) THEN
C- Compute Qsat and d.Qsat/d.T :
DO J=1,NGP
QSAT(J)=0.
sigP = PS(1)
IF (SIG.GT.0.0) sigP=SIG*PS(J)
IF (TA(J).GE.T0) THEN
tmpQ = E0*EXP(C1*(TA(J)-T0)/(TA(J)-T1))
QSAT(J)= QS1*tmpQ/(sigP-QS2*tmpQ)
recT = 1. _d 0 / (TA(J)-T1)
RH(J) = QSAT(J)*C1*(T0-T1)*recT*recT*sigP/(sigP-QS2*tmpQ)
ELSE IF ( TA(J).GT.T2) THEN
tmpQ = E0*EXP(C2*(TA(J)-T0)/(TA(J)-T2))
QSAT(J)= QS1*tmpQ/(sigP-QS2*tmpQ)
recT = 1. _d 0 / (TA(J)-T2)
RH(J) = QSAT(J)*C2*(T0-T2)*recT*recT*sigP/(sigP-QS2*tmpQ)
ENDIF
ENDDO
RETURN
ENDIF
DO 110 J=1,NGP
QSAT(J)=0.
IF (TA(J).GE.T0) THEN
QSAT(J)=E0*EXP(C1*(TA(J)-T0)/(TA(J)-T1))
ELSE IF ( TA(J).GT.T2) THEN
QSAT(J)=E0*EXP(C2*(TA(J)-T0)/(TA(J)-T2))
ENDIF
110 CONTINUE
C
IF (SIG.LE.0.0) THEN
DO 120 J=1,NGP
QSAT(J)= QS1*QSAT(J)/( PS(1) - QS2*QSAT(J))
120 CONTINUE
ELSE
DO 130 J=1,NGP
QSAT(J)= QS1*QSAT(J)/(SIG*PS(J)- QS2*QSAT(J))
130 CONTINUE
ENDIF
C
C--- 2. Compute rel.hum. RH=Q/Qsat (IMODE>0), or Q=RH*Qsat (IMODE<0)
C
IF (IMODE.GT.0) THEN
DO 210 J=1,NGP
IF(QSAT(J).NE.0.) then
RH(J)=QA(J)/QSAT(J)
ELSE
RH(J)=0.
ENDIF
210 CONTINUE
ELSE IF (IMODE.LT.0) THEN
DO 220 J=1,NGP
QA(J)=RH(J)*QSAT(J)
220 CONTINUE
ENDIF
#endif /* ALLOW_AIM */
RETURN
END
SUBROUTINE ZMEDDY (NLON,NLAT,FF,ZM,EDDY)
IMPLICIT NONE
C *** Decompose a field into zonal-mean and eddy component
C-- Routine arguments:
INTEGER NLON, NLAT
_RL FF(NLON,NLAT), ZM(NLAT), EDDY(NLON,NLAT)
#ifdef ALLOW_AIM
C-- Local variables:
INTEGER I,J
C- jmc: declare all local variables:
_RL RNLON
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
RNLON=1./NLON
C
DO 130 J=1,NLAT
C
ZM(J)=0.
DO 110 I=1,NLON
ZM(J)=ZM(J)+FF(I,J)
110 CONTINUE
ZM(J)=ZM(J)*RNLON
C
DO 120 I=1,NLON
EDDY(I,J)=FF(I,J)-ZM(J)
120 CONTINUE
C
130 CONTINUE
C
C--
#endif /* ALLOW_AIM */
RETURN
END