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