1 |
C $Header: /u/gcmpack/models/MITgcmUV/pkg/aim/phy_suflux.F,v 1.3 2001/05/29 19:28:53 cnh Exp $ |
2 |
C $Name: $ |
3 |
|
4 |
SUBROUTINE SUFLUX (PSA,UA,VA,TA,QA,RH,QSAT,Vsurfsq,PHI, |
5 |
* PHI0,FMASK,TLAND,TSEA,SWAV,SSR,SLR, |
6 |
* DRAG,USTR,VSTR,SHF,EVAP,T0,Q0,QSAT0,SPEED0, |
7 |
& myThid) |
8 |
C-- |
9 |
C-- SUBROUTINE SUFLUX (PSA,UA,VA,TA,QA,RH,PHI, |
10 |
C-- * PHI0,FMASK,TLAND,TSEA,SWAV, |
11 |
C-- * USTR,VSTR,SHF,EVAP) |
12 |
C-- |
13 |
C-- Purpose: Compute surface fluxes of momentum, energy and moisture |
14 |
C-- Input: PSA = norm. surface pressure [p/p0] (2-dim) |
15 |
C-- UA = u-wind (3-dim) |
16 |
C-- VA = v-wind (3-dim) |
17 |
C-- TA = temperature (3-dim) |
18 |
C-- QA = specific humidity [g/kg] (3-dim) |
19 |
C-- RH = relative humidity (3-dim) |
20 |
C-- PHI = geopotential (3-dim) |
21 |
C-- PHI0 = surface geopotential (2-dim) |
22 |
C-- FMASK = fractional land-sea mask (2-dim) |
23 |
C-- TLAND = land-surface temperature (2-dim) |
24 |
C-- TSEA = sea-surface temperature (2-dim) |
25 |
C-- SWET = soil wetness availability [0-1] (2-dim) |
26 |
C-- Output: USTR = u stress (2-dim) |
27 |
C-- VSTR = v stress (2-dim) |
28 |
C-- SHF = sensible heat flux (2-dim) |
29 |
C-- EVAP = evaporation [g/(m^2 s)] (2-dim) |
30 |
C - added (jmc) : |
31 |
C-- Vsurfsq = square of surface wind speed (2-dim,input) |
32 |
C-- ==> UA,VA are no longer used |
33 |
C-- DRAG = surface Drag term (= Cd*|V|) (2-dim,output) |
34 |
C-- ==> USTR,VSTR are no longer used |
35 |
C-- myThid = Instance number of this instance of the |
36 |
C-- the routine. |
37 |
C-- |
38 |
|
39 |
|
40 |
IMPLICIT rEAL*8 (A-H,O-Z) |
41 |
|
42 |
|
43 |
C Resolution parameters |
44 |
#include "EEPARAMS.h" |
45 |
#include "atparam.h" |
46 |
#include "atparam1.h" |
47 |
#include "Lev_def.h" |
48 |
C |
49 |
INTEGER NLON,NLAT,NLEV,NGP |
50 |
PARAMETER ( NLON=IX, NLAT=IL, NLEV=KX, NGP=NLON*NLAT ) |
51 |
|
52 |
C Physical constants + functions of sigma and latitude |
53 |
|
54 |
#include "com_physcon.h" |
55 |
|
56 |
C Surface flux constants |
57 |
|
58 |
#include "com_sflcon.h" |
59 |
|
60 |
C |
61 |
REAL PSA(NGP), UA(NGP,NLEV), VA(NGP,NLEV), TA(NGP,NLEV), |
62 |
* QA(NGP,NLEV), RH(NGP,NLEV), QSAT(NGP,NLEV), PHI(NGP,NLEV), |
63 |
* PHI0(NGP), FMASK(NGP), TLAND(NGP), TSEA(NGP), SWAV(NGP), |
64 |
* SSR(NGP), SLR(NGP) |
65 |
REAL Vsurfsq(NGP), DRAG(NGP) |
66 |
INTEGER myThid |
67 |
C |
68 |
REAL USTR(NGP,3), VSTR(NGP,3), SHF(NGP,3), EVAP(NGP,3) |
69 |
C |
70 |
REAL U0(NGP),V0(NGP),T0(NGP,2),Q0(NGP),QSAT0(NGP,2),DENVV(NGP,3) |
71 |
REAL SPEED0(NGP) |
72 |
REAL WORK(NGP) |
73 |
INTEGER NL1(NGP) |
74 |
REAL AUX(NGP) |
75 |
C |
76 |
REAL pSurfs(NLEV) |
77 |
DATA pSurfs /75 _d 2,250 _d 2,500 _d 2,775 _d 2,950 _d 2 / |
78 |
C |
79 |
REAL pSurfw(NLEV) |
80 |
DATA pSurfw /150 _d 2,350 _d 2,650 _d 2,900 _d 2,1000 _d 2/ |
81 |
Cchdbg |
82 |
INTEGER NPAS |
83 |
SAVE NPAS |
84 |
LOGICAL Ifirst |
85 |
SAVE Ifirst |
86 |
DATA Ifirst /.TRUE./ |
87 |
REAL T0moy1(NGP) |
88 |
REAL T0moy2(NGP) |
89 |
SAVE T0moy2, T0moy1 |
90 |
REAL denmoy1(NGP) |
91 |
REAL denmoy2(NGP) |
92 |
SAVE denmoy1, denmoy2 |
93 |
REAL Q0moy(NGP) |
94 |
SAVE Q0moy |
95 |
REAL factwind2 |
96 |
INTEGER J |
97 |
C |
98 |
C-- 1. Extrapolation of wind, temp, hum. and density to the surface |
99 |
C |
100 |
C 1.1 Wind components |
101 |
C |
102 |
DO J=1,NGP |
103 |
c_jmc IF ( NLEVxyU(J) .GT. 0 ) THEN |
104 |
c_jmc U0(J) = FWIND0*UA(J,NLEVxyU(J)) |
105 |
c_jmc ENDIF |
106 |
c_jmc IF ( NLEVxyV(J) .GT. 0 ) THEN |
107 |
c_jmc V0(J) = FWIND0*VA(J,NLEVxyV(J)) |
108 |
c_jmc ENDIF |
109 |
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
110 |
C-- Note(jmc): If we interpolate U,V from C_grid to A_grid (center of the mesh) |
111 |
C the above distinction between NLEVxy_Upt,Vpt is no longer valid ! |
112 |
C same thing (again) few lines below. |
113 |
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
114 |
IF ( NLEVxy(J,myThid) .GT. 0 ) THEN |
115 |
U0(J) = FWIND0*UA(J,NLEVxy(J,myThid)) |
116 |
V0(J) = FWIND0*VA(J,NLEVxy(J,myThid)) |
117 |
ELSE |
118 |
U0(J) = 0.0 |
119 |
V0(J) = 0.0 |
120 |
ENDIF |
121 |
C U0(J) = UA(J,5) |
122 |
C V0(J) = VA(J,5) |
123 |
ENDDO |
124 |
C |
125 |
C 1.2 Temperature |
126 |
C |
127 |
GTEMP0 = 1.D0-FTEMP0 |
128 |
RCP = 1.D0/CP |
129 |
DO J=1,NGP |
130 |
NL1(J)=NLEVxy(J,myThid)-1 |
131 |
ENDDO |
132 |
C |
133 |
DO J=1,NGP |
134 |
IF ( NLEVxy(J,myThid) .GT. 0 ) THEN |
135 |
T0(J,1) = TA(J,NLEVxy(J,myThid))+WVI(NLEVxy(J,myThid),2)* |
136 |
& (TA(J,NLEVxy(J,myThid))-TA(J,NL1(J))) |
137 |
ccccc T0(J,2) = TA(J,NLEVxy(J))+RCP*(PHI(J,NLEVxy(J))-PHI0(J)) |
138 |
T0(J,2) = TA(J,NLEVxy(J,myThid))* |
139 |
& ((Psurfw(NLEVxy(J,myThid))/ |
140 |
& Psurfs(Nlevxy(J,myThid)))**(RD/CP)) |
141 |
ELSE |
142 |
T0(J,1) = 273.16 _d 0 |
143 |
T0(J,2) = 273.16 _d 0 |
144 |
ENDIF |
145 |
ENDDO |
146 |
C |
147 |
DO J=1,NGP |
148 |
T0(J,1) = FTEMP0*T0(J,1)+GTEMP0*T0(J,2) |
149 |
ENDDO |
150 |
C |
151 |
C 1.3 Spec. humidity |
152 |
C |
153 |
GHUM0 = 1. _d 0-FHUM0 |
154 |
C |
155 |
DO J=1,NGP |
156 |
IF ( NLEVxy(J,myThid) .GT. 0 ) THEN |
157 |
WORK(J)=RH(J,Nlevxy(J,myThid)) |
158 |
cchdbg |
159 |
c WORK(J) = RH(J,NLEVxy(J))+WVI(NLEVxy(J),2)* |
160 |
c & (RH(J,NLEVxy(J))-RH(J,NL1(J))) |
161 |
cchdbg |
162 |
ENDIF |
163 |
ENDDO |
164 |
C |
165 |
CALL SHTORH (-1,NGP,T0,PSA,1. _d 0,Q0,WORK,QSAT0,myThid) |
166 |
C |
167 |
DO J=1,NGP |
168 |
IF ( NLEVxy(J,myThid) .GT. 0 ) THEN |
169 |
Q0(J)=FHUM0*Q0(J)+GHUM0*QA(J,NLEVxy(J,myThid)) |
170 |
ENDIF |
171 |
ENDDO |
172 |
|
173 |
C 1.4 Density * wind speed (including gustiness factor) |
174 |
|
175 |
PRD = P0/RD |
176 |
VG2 = VGUST*VGUST |
177 |
factwind2 = FWIND0*FWIND0 |
178 |
C |
179 |
DO J=1,NGP |
180 |
c_jmc SPEED0(J)=SQRT(U0(J)*U0(J)+V0(J)*V0(J)+VG2) |
181 |
SPEED0(J)=SQRT(factwind2*Vsurfsq(J)+VG2) |
182 |
DENVV(J,3)=(PRD*PSA(J)/T0(J,1))*SPEED0(J) |
183 |
c_jmc DENVV(J,3)=(PRD*PSA(J)/T0(J,1))* |
184 |
c_jmc* SQRT(U0(J)*U0(J)+V0(J)*V0(J)+VG2) |
185 |
cchdbg DENVV(J,3)=0.7*(PRD*PSA(J)/T0(J,1))* |
186 |
cchdbg * SQRT(U0(J)*U0(J)+V0(J)*V0(J)+VG2) |
187 |
ENDDO |
188 |
C |
189 |
C 1.5 Stability correction for heat/moisture fluxes |
190 |
C |
191 |
DTHETAF = 3. _d 0 |
192 |
FSTAB = 0.67 _d 0 |
193 |
RDTH = FSTAB/DTHETAF |
194 |
C |
195 |
DO J=1,NGP |
196 |
FSLAND=1.+MIN(DTHETAF,MAX(-DTHETAF,TLAND(J)-T0(J,2)))*RDTH |
197 |
FSSEA =1.+MIN(DTHETAF,MAX(-DTHETAF, TSEA(J)-T0(J,2)))*RDTH |
198 |
aux(J)=FSLAND |
199 |
DENVV(J,1)=DENVV(J,3)*FSLAND |
200 |
DENVV(J,2)=DENVV(J,3)*FSSEA |
201 |
cchdbg DENVV(J,1)=DENVV(J,3) |
202 |
cchdbg DENVV(J,2)=DENVV(J,3) |
203 |
ENDDO |
204 |
C |
205 |
C-- 2. Computation of fluxes over land and sea |
206 |
C |
207 |
C 2.1 Wind stress |
208 |
C |
209 |
DO J=1,NGP |
210 |
c_jmc IF ( NLEVxyu(J) .GT. 0 ) THEN |
211 |
c_jmc USTR(J,1) = -CDL*DENVV(J,3)*UA(J,NLEVxyu(J)) |
212 |
c_jmc USTR(J,2) = -CDS*DENVV(J,3)*UA(J,NLEVxyu(J)) |
213 |
c_jmc ENDIF |
214 |
c_jmc IF ( NLEVxyv(J) .GT. 0 ) THEN |
215 |
c_jmc VSTR(J,1) = -CDL*DENVV(J,3)*VA(J,NLEVxyv(J)) |
216 |
c_jmc VSTR(J,2) = -CDS*DENVV(J,3)*VA(J,NLEVxyv(j)) |
217 |
c_jmc ENDIF |
218 |
C-- cf Note(jmc) above. |
219 |
IF ( NLEVxy(J,myThid) .GT. 0 ) THEN |
220 |
USTR(J,1) = -CDL*DENVV(J,3)*UA(J,NLEVxy(J,myThid)) |
221 |
USTR(J,2) = -CDS*DENVV(J,3)*UA(J,NLEVxy(J,myThid)) |
222 |
VSTR(J,1) = -CDL*DENVV(J,3)*VA(J,NLEVxy(J,myThid)) |
223 |
VSTR(J,2) = -CDS*DENVV(J,3)*VA(J,NLEVxy(J,myThid)) |
224 |
ELSE |
225 |
USTR(J,1) = 0. |
226 |
USTR(J,2) = 0. |
227 |
VSTR(J,1) = 0. |
228 |
VSTR(J,2) = 0. |
229 |
ENDIF |
230 |
ENDDO |
231 |
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
232 |
c Compute surface drag term (= C_drag*|V| ) to allow direct computation |
233 |
c of surface stress on C-grid. |
234 |
c add Land + Sea contributions ; Convert to surface pressure level |
235 |
DO J=1,NGP |
236 |
IF ( NLEVxy(J,myThid) .GT. 0 ) THEN |
237 |
DRAG(J) = ( CDS+FMASK(J)*(CDL-CDS) ) |
238 |
& * DENVV(J,3) * GRDSIG(NLEVxy(J,myThid)) |
239 |
ELSE |
240 |
DRAG(J) = 0. |
241 |
ENDIF |
242 |
ENDDO |
243 |
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
244 |
C |
245 |
C 2.2 Sensible heat flux (from clim. TS over land) |
246 |
C |
247 |
CHLCP = CHL*CP |
248 |
CHSCP = CHS*CP |
249 |
C |
250 |
DO J=1,NGP |
251 |
SHF(J,1) = CHLCP*DENVV(J,1)*(TLAND(J)-T0(J,1)) |
252 |
SHF(J,2) = CHSCP*DENVV(J,2)*(TSEA(J) -T0(J,1)) |
253 |
ENDDO |
254 |
C |
255 |
C **************************************************** |
256 |
cchdbg |
257 |
IF (Ifirst) then |
258 |
npas=0. |
259 |
do J=1,NGP |
260 |
T0moy1(J)=0. |
261 |
T0moy2(J)=0. |
262 |
denmoy1(J)=0. |
263 |
denmoy2(J)=0. |
264 |
Q0moy(J)=0. |
265 |
enddo |
266 |
ifirst=.false. |
267 |
ENDIF |
268 |
C |
269 |
npas=npas+1 |
270 |
DO J=1,NGP |
271 |
T0moy1(J)=T0moy1(J) + T0(J,1) |
272 |
T0moy2(J)=T0moy2(J) + T0(J,2) |
273 |
denmoy1(J)=denmoy1(J) + DENVV(J,1) |
274 |
denmoy2(J)=denmoy2(J) + DENVV(J,2) |
275 |
Q0moy(J)=Q0moy(J) + Q0(J) |
276 |
ENDDO |
277 |
C |
278 |
if(npas.eq.5760) then |
279 |
DO J=1,NGP |
280 |
T0moy1(J)=T0moy1(J)/float(npas) |
281 |
T0moy2(J)=T0moy2(J)/float(npas) |
282 |
denmoy1(J)=denmoy1(J)/float(npas) |
283 |
denmoy2(J)=denmoy2(J)/float(npas) |
284 |
Q0moy(J)=Q0moy(J)/float(npas) |
285 |
ENDDO |
286 |
C |
287 |
open(73,file='Tmoy1',form='unformatted') |
288 |
write(73) T0moy1 |
289 |
close(73) |
290 |
open(74,file='Tmoy2',form='unformatted') |
291 |
write(74) T0moy2 |
292 |
close(74) |
293 |
open(73,file='denmoy1',form='unformatted') |
294 |
write(73) denmoy1 |
295 |
close(73) |
296 |
open(74,file='denmoy2',form='unformatted') |
297 |
write(74) denmoy2 |
298 |
close(74) |
299 |
open(74,file='Q0moy',form='unformatted') |
300 |
write(74) Q0moy |
301 |
close(74) |
302 |
ENDIF |
303 |
C |
304 |
cchdbg |
305 |
C **************************************************** |
306 |
C |
307 |
c CALL DUMP_WRITE2D ( NGP,1,'T0.', nIter,T0 , iErr) |
308 |
c CALL DUMP_WRITE2D ( NGP,3,'DENVV.', nIter,DENVV , iErr) |
309 |
c CALL DUMP_WRITE2D ( NGP,1,'TSEA.', nIter,TSEA , iErr) |
310 |
c CALL DUMP_WRITE2D ( NGP,1,'TLAND.', nIter,TLAND , iErr) |
311 |
c CALL DUMP_WRITE2D ( NGP,1,'AUX.', nIter,aux , iErr) |
312 |
C |
313 |
C 2.3 Evaporation |
314 |
C |
315 |
CALL SHTORH (0,NGP,TLAND,PSA,1. _d 0,QDUMMY,RDUMMY,QSAT0(1,1),myThid) |
316 |
CALL SHTORH (0,NGP,TSEA ,PSA,1. _d 0,QDUMMY,RDUMMY,QSAT0(1,2),myThid) |
317 |
|
318 |
DO J=1,NGP |
319 |
Cdj EVAP(J,1) = CHL*DENVV(J,1)*MAX(0. _d 0,SWAV(J)*QSAT0(J,1)-Q0(J)) |
320 |
c_jmc: EVAP(J,1)=CHL*DENVV(J,1)*MIN(1. _d 0,SWAV(J))*(QSAT0(J,1)-Q0(J)) |
321 |
EVAP(J,1)=CHL*DENVV(J,1)*MAX(0. _d 0,SWAV(J)*(QSAT0(J,1)-Q0(J))) |
322 |
cdj try new scheme : assume that the net heat flux on land is zero |
323 |
cdj at all time and at all points (this is equivalent |
324 |
cdj to say that the land has a zero heat capacity) |
325 |
cdj EVAP(J,1) = SSR(J)-SLR(J)-SHF(J,1) |
326 |
c_jmc: EVAP(J,2) = CHS*DENVV(J,2)*(QSAT0(J,2)-Q0(J)) |
327 |
EVAP(J,2) = CHS*DENVV(J,2)*MAX(0. _d 0,QSAT0(J,2)-Q0(J)) |
328 |
ENDDO |
329 |
|
330 |
|
331 |
C-- 3. Weighted average of fluxes according to land-sea mask |
332 |
|
333 |
DO J=1,NGP |
334 |
USTR(J,3) = USTR(J,2)+FMASK(J)*(USTR(J,1)-USTR(J,2)) |
335 |
VSTR(J,3) = VSTR(J,2)+FMASK(J)*(VSTR(J,1)-VSTR(J,2)) |
336 |
SHF(J,3) = SHF(J,2)+FMASK(J)*( SHF(J,1)- SHF(J,2)) |
337 |
EVAP(J,3) = EVAP(J,2)+FMASK(J)*(EVAP(J,1)-EVAP(J,2)) |
338 |
cdj |
339 |
QSAT0(J,1) = QSAT0(J,2)+FMASK(J)*(QSAT0(J,1)-QSAT0(J,2)) |
340 |
cdj |
341 |
ENDDO |
342 |
|
343 |
C-- |
344 |
RETURN |
345 |
END |