/[MITgcm]/MITgcm/pkg/aim/phy_suflux.F
ViewVC logotype

Annotation of /MITgcm/pkg/aim/phy_suflux.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.3 - (hide annotations) (download)
Tue May 29 19:28:53 2001 UTC (23 years ago) by cnh
Branch: MAIN
CVS Tags: checkpoint40pre1, checkpoint40pre2, checkpoint40pre5, checkpoint40pre6, checkpoint40pre8, checkpoint40pre4, checkpoint40pre3, checkpoint40pre7
Changes since 1.2: +89 -33 lines
Updates for multi-threaded AIM with support for both latlon
and CS.
Needs compatible changes to verfication/

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

  ViewVC Help
Powered by ViewVC 1.1.22