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

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

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

revision 1.1 by cnh, Fri Jan 26 00:14:32 2001 UTC revision 1.2 by adcroft, Fri Feb 2 21:36:29 2001 UTC
# Line 0  Line 1 
1    C $Header$
2    C $Name$
3    
4          SUBROUTINE SUFLUX (PSA,UA,VA,TA,QA,RH,QSAT,PHI,
5         *                   PHI0,FMASK,TLAND,TSEA,SWAV,SSR,SLR,
6         *                   USTR,VSTR,SHF,EVAP,T0,Q0,QSAT0,SPEED0)
7    C--
8    C--   SUBROUTINE SUFLUX (PSA,UA,VA,TA,QA,RH,PHI,
9    C--  *                   PHI0,FMASK,TLAND,TSEA,SWAV,
10    C--  *                   USTR,VSTR,SHF,EVAP)
11    C--
12    C--   Purpose: Compute surface fluxes of momentum, energy and moisture
13    C--   Input:   PSA    = norm. surface pressure [p/p0] (2-dim)
14    C--            UA     = u-wind                        (3-dim)
15    C--            VA     = v-wind                        (3-dim)
16    C--            TA     = temperature                   (3-dim)
17    C--            QA     = specific humidity [g/kg]      (3-dim)
18    C--            RH     = relative humidity             (3-dim)
19    C--            PHI    = geopotential                  (3-dim)
20    C--            PHI0   = surface geopotential          (2-dim)
21    C--            FMASK  = fractional land-sea mask         (2-dim)
22    C--            TLAND  = land-surface temperature         (2-dim)
23    C--            TSEA   =  sea-surface temperature         (2-dim)
24    C--            SWET   = soil wetness availability [0-1]  (2-dim)
25    C--   Output:  USTR   = u stress                         (2-dim)
26    C--            VSTR   = v stress                         (2-dim)
27    C--            SHF    = sensible heat flux               (2-dim)
28    C--            EVAP   = evaporation [g/(m^2 s)]          (2-dim)
29    C--
30    
31    
32          IMPLICIT rEAL*8 (A-H,O-Z)
33    
34    
35    C     Resolution parameters
36    #include "atparam.h"
37    #include "atparam1.h"
38    #include "Lev_def.h"
39    C
40          PARAMETER ( NLON=IX, NLAT=IL, NLEV=KX, NGP=NLON*NLAT )
41    
42    C     Physical constants + functions of sigma and latitude
43    
44    #include "com_physcon.h"
45    
46    C     Surface flux constants
47    
48    #include "com_sflcon.h"
49    C
50          REAL PSA(NGP), UA(NGP,NLEV), VA(NGP,NLEV), TA(NGP,NLEV),
51         *     QA(NGP,NLEV), RH(NGP,NLEV), QSAT(NGP,NLEV), PHI(NGP,NLEV),
52         *     PHI0(NGP), FMASK(NGP), TLAND(NGP), TSEA(NGP), SWAV(NGP),
53         *     SSR(NGP), SLR(NGP)
54    C
55          REAL USTR(NGP,3), VSTR(NGP,3), SHF(NGP,3), EVAP(NGP,3)
56    C                                                                      
57          REAL U0(NGP),V0(NGP),T0(NGP,2),Q0(NGP),QSAT0(NGP,2),DENVV(NGP,3)
58          REAL SPEED0(NGP)
59          REAL WORK(NGP)
60          INTEGER NL1(NGP)
61          REAL AUX(NGP)
62    C
63          REAL    pSurfs(NLEV)
64          DATA    pSurfs /75 _d 2,250 _d 2,500 _d 2,775 _d 2,950 _d 2 /
65    C
66          REAL    pSurfw(NLEV)
67          DATA    pSurfw /150 _d 2,350 _d 2,650 _d 2,900 _d 2,1000 _d 2/
68    Cchdbg
69          INTEGER NPAS
70          SAVE NPAS
71          LOGICAL Ifirst
72          SAVE Ifirst
73          DATA Ifirst /.TRUE./
74          REAL T0moy1(NGP)
75          REAL T0moy2(NGP)
76          SAVE T0moy2, T0moy1
77          REAL denmoy1(NGP)
78          REAL denmoy2(NGP)
79          SAVE denmoy1, denmoy2
80          REAL Q0moy(NGP)
81          SAVE Q0moy
82    C
83    C--   1. Extrapolation of wind, temp, hum. and density to the surface
84    C
85    C     1.1 Wind components
86    C  
87          DO J=1,NGP
88            IF ( NLEVxyU(J) .GT. 0 ) THEN
89             U0(J) = FWIND0*UA(J,NLEVxyU(J))
90            ENDIF
91            IF ( NLEVxyV(J) .GT. 0 ) THEN
92             V0(J) = FWIND0*VA(J,NLEVxyV(J))
93            ENDIF
94    C       U0(J) = UA(J,5)
95    C       V0(J) = VA(J,5)
96          ENDDO
97    C
98    C     1.2 Temperature
99    C
100          GTEMP0 = 1.D0-FTEMP0
101          RCP = 1.D0/CP
102          DO J=1,NGP
103            NL1(J)=NLEVxy(J)-1
104          ENDDO
105    C
106          DO J=1,NGP
107           IF ( NLEVxy(J) .GT. 0 ) THEN
108            T0(J,1) = TA(J,NLEVxy(J))+WVI(NLEVxy(J),2)*
109         &                          (TA(J,NLEVxy(J))-TA(J,NL1(J)))
110    ccccc        T0(J,2) = TA(J,NLEVxy(J))+RCP*(PHI(J,NLEVxy(J))-PHI0(J))
111            T0(J,2) = TA(J,NLEVxy(J))*
112         &         ((Psurfw(NLEVxy(J))/Psurfs(Nlevxy(J)))**(RD/CP))
113           ELSE
114            T0(J,1) = 273.16 _d 0
115            T0(J,2) = 273.16 _d 0
116           ENDIF
117          ENDDO
118    C
119          DO J=1,NGP
120            T0(J,1) = FTEMP0*T0(J,1)+GTEMP0*T0(J,2)
121          ENDDO
122    C
123    C     1.3 Spec. humidity
124    C
125          GHUM0 = 1. _d 0-FHUM0
126    C
127          DO J=1,NGP
128           IF ( NLEVxy(J) .GT. 0 ) THEN
129            WORK(J)=RH(J,Nlevxy(J))
130    cchdbg
131    c        WORK(J) = RH(J,NLEVxy(J))+WVI(NLEVxy(J),2)*
132    c     &                          (RH(J,NLEVxy(J))-RH(J,NL1(J)))
133    cchdbg
134           ENDIF
135          ENDDO
136    C
137          CALL SHTORH (-1,NGP,T0,PSA,1. _d 0,Q0,WORK,QSAT0)
138    C
139          DO J=1,NGP
140           IF ( NLEVxy(J) .GT. 0 ) THEN
141            Q0(J)=FHUM0*Q0(J)+GHUM0*QA(J,NLEVxy(J))
142           ENDIF
143          ENDDO
144    
145    C     1.4 Density * wind speed (including gustiness factor)
146    
147          PRD = P0/RD
148          VG2 = VGUST*VGUST
149    C
150          DO J=1,NGP
151            DENVV(J,3)=(PRD*PSA(J)/T0(J,1))*
152         *           SQRT(U0(J)*U0(J)+V0(J)*V0(J)+VG2)
153    cchdbg        DENVV(J,3)=0.7*(PRD*PSA(J)/T0(J,1))*
154    cchdbg     *           SQRT(U0(J)*U0(J)+V0(J)*V0(J)+VG2)
155            SPEED0(J)=SQRT(U0(J)*U0(J)+V0(J)*V0(J)+VG2)
156          ENDDO
157    C
158    C     1.5 Stability correction for heat/moisture fluxes
159    C
160          DTHETAF = 3. _d 0
161          FSTAB = 0.67 _d 0
162          RDTH = FSTAB/DTHETAF
163    C
164          DO J=1,NGP
165            FSLAND=1.+MIN(DTHETAF,MAX(-DTHETAF,TLAND(J)-T0(J,2)))*RDTH
166            FSSEA =1.+MIN(DTHETAF,MAX(-DTHETAF, TSEA(J)-T0(J,2)))*RDTH
167            aux(J)=FSLAND
168            DENVV(J,1)=DENVV(J,3)*FSLAND
169            DENVV(J,2)=DENVV(J,3)*FSSEA
170    cchdbg        DENVV(J,1)=DENVV(J,3)
171    cchdbg        DENVV(J,2)=DENVV(J,3)
172          ENDDO
173    C
174    C--   2. Computation of fluxes over land and sea
175    C
176    C     2.1 Wind stress
177    C
178          DO J=1,NGP
179           IF ( NLEVxyu(J) .GT. 0  ) THEN
180            USTR(J,1) = -CDL*DENVV(J,3)*UA(J,NLEVxyu(J))
181            USTR(J,2) = -CDS*DENVV(J,3)*UA(J,NLEVxyu(J))
182           ENDIF
183           IF ( NLEVxyv(J) .GT. 0  ) THEN
184            VSTR(J,1) = -CDL*DENVV(J,3)*VA(J,NLEVxyv(J))
185            VSTR(J,2) = -CDS*DENVV(J,3)*VA(J,NLEVxyv(j))
186           ENDIF
187          ENDDO
188    C
189    C     2.2 Sensible heat flux (from clim. TS over land)
190    C
191          CHLCP = CHL*CP
192          CHSCP = CHS*CP
193    C
194          DO J=1,NGP
195            SHF(J,1) = CHLCP*DENVV(J,1)*(TLAND(J)-T0(J,1))
196            SHF(J,2) = CHSCP*DENVV(J,2)*(TSEA(J) -T0(J,1))
197          ENDDO
198    C
199    C ****************************************************
200    cchdbg
201          IF (Ifirst) then
202            npas=0.
203            do J=1,NGP
204              T0moy1(J)=0.
205              T0moy2(J)=0.
206              denmoy1(J)=0.
207              denmoy2(J)=0.
208              Q0moy(J)=0.
209            enddo
210            ifirst=.false.
211          ENDIF
212    C
213          npas=npas+1
214          DO J=1,NGP
215            T0moy1(J)=T0moy1(J) + T0(J,1)
216            T0moy2(J)=T0moy2(J) + T0(J,2)
217            denmoy1(J)=denmoy1(J) + DENVV(J,1)
218            denmoy2(J)=denmoy2(J) + DENVV(J,2)
219            Q0moy(J)=Q0moy(J) + Q0(J)
220          ENDDO
221    C
222          if(npas.eq.5760) then
223            DO J=1,NGP
224              T0moy1(J)=T0moy1(J)/float(npas)
225              T0moy2(J)=T0moy2(J)/float(npas)
226              denmoy1(J)=denmoy1(J)/float(npas)
227              denmoy2(J)=denmoy2(J)/float(npas)
228              Q0moy(J)=Q0moy(J)/float(npas)
229            ENDDO
230    C
231            open(73,file='Tmoy1',form='unformatted')
232            write(73) T0moy1
233            close(73)
234            open(74,file='Tmoy2',form='unformatted')
235            write(74) T0moy2
236            close(74)
237            open(73,file='denmoy1',form='unformatted')
238            write(73) denmoy1
239            close(73)
240            open(74,file='denmoy2',form='unformatted')
241            write(74) denmoy2
242            close(74)
243            open(74,file='Q0moy',form='unformatted')
244            write(74) Q0moy
245            close(74)
246          ENDIF
247    C
248    cchdbg
249    C ****************************************************
250    C
251    c        CALL DUMP_WRITE2D ( NGP,1,'T0.', nIter,T0 , iErr)
252    c        CALL DUMP_WRITE2D ( NGP,3,'DENVV.', nIter,DENVV , iErr)
253    c        CALL DUMP_WRITE2D ( NGP,1,'TSEA.', nIter,TSEA , iErr)
254    c        CALL DUMP_WRITE2D ( NGP,1,'TLAND.', nIter,TLAND , iErr)
255    c        CALL DUMP_WRITE2D ( NGP,1,'AUX.', nIter,aux , iErr)
256    C
257    C     2.3 Evaporation
258    C
259          CALL SHTORH (0,NGP,TLAND,PSA,1. _d 0,QDUMMY,RDUMMY,QSAT0(1,1))
260          CALL SHTORH (0,NGP,TSEA ,PSA,1. _d 0,QDUMMY,RDUMMY,QSAT0(1,2))
261    
262          DO J=1,NGP
263    Cdj        EVAP(J,1) = CHL*DENVV(J,1)*MAX(0. _d 0,SWAV(J)*QSAT0(J,1)-Q0(J))
264            EVAP(J,1)=CHL*DENVV(J,1)*MIN(1. _d 0,SWAV(J))*(QSAT0(J,1)-Q0(J))
265    cdj try new scheme : assume that the net heat flux on land is zero
266    cdj                  at all time and at all points (this is equivalent
267    cdj                  to say that the land has a zero heat capacity)
268    cdj        EVAP(J,1) = SSR(J)-SLR(J)-SHF(J,1)
269            EVAP(J,2) = CHS*DENVV(J,2)*(QSAT0(J,2)-Q0(J))
270          ENDDO
271    
272    
273    C--   3. Weighted average of fluxes according to land-sea mask
274    
275          DO J=1,NGP
276            USTR(J,3) = USTR(J,2)+FMASK(J)*(USTR(J,1)-USTR(J,2))
277            VSTR(J,3) = VSTR(J,2)+FMASK(J)*(VSTR(J,1)-VSTR(J,2))
278             SHF(J,3) =  SHF(J,2)+FMASK(J)*( SHF(J,1)- SHF(J,2))
279            EVAP(J,3) = EVAP(J,2)+FMASK(J)*(EVAP(J,1)-EVAP(J,2))
280    cdj
281            QSAT0(J,1) = QSAT0(J,2)+FMASK(J)*(QSAT0(J,1)-QSAT0(J,2))
282    cdj
283          ENDDO
284    
285    C--
286          RETURN
287          END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22