/[MITgcm]/MITgcm_contrib/jscott/igsm/src/sur4clm.F
ViewVC logotype

Contents of /MITgcm_contrib/jscott/igsm/src/sur4clm.F

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


Revision 1.2 - (show annotations) (download)
Mon Apr 23 21:20:18 2007 UTC (18 years, 3 months ago) by jscott
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +28 -19 lines
bring igsm atmos code up to date

1
2 #include "ctrparam.h"
3
4 ! ==========================================================
5 !
6 ! SURFACE.F: THIS SUBROUTINE CALCULATES THE SURFACE FLUXES
7 ! WHICH INCLUDE SENSIBLE HEAT, EVAPORATION,
8 ! THERMAL RADIATION, AND MOMENTUM DRAG. IT ALSO
9 ! CALCULATES INSTANTANEOUSLY SURFACE TEMPERATURE,
10 ! SURFACE SPECIFIC HUMIDITY, AND SURFACE WIND
11 ! COMPONENTS.
12 !
13 ! ----------------------------------------------------------
14 !
15 ! Author of Chemistry Modules: Chien Wang
16 !
17 ! ----------------------------------------------------------
18 !
19 ! Revision History:
20 !
21 ! When Who What
22 ! ---- ---------- -------
23 ! 073100 Chien Wang repack based on CliChem3 and add cpp
24 ! 092301 Chien Wang add bc and oc
25 !
26 ! ==========================================================
27
28 SUBROUTINE SUR4CLM
29
30 C**** 5802.
31 C**** THIS SUBROUTINE CALCULATES THE SURFACE FLUXES WHICH INCLUDE 5803.
32 C**** SENSIBLE HEAT, EVAPORATION, THERMAL RADIATION, AND MOMENTUM 5804.
33 C**** DRAG. IT ALSO CALCULATES INSTANTANEOUSLY SURFACE TEMPERATURE, 5805.
34 C**** SURFACE SPECIFIC HUMIDITY, AND SURFACE WIND COMPONENTS. 5806.
35 C**** 5807.
36
37 #if ( defined CLM )
38 #if ( defined CPL_CHEM )
39 !
40 #include "chem_para"
41 #include "chem_com"
42 !
43 #endif
44
45 #include "BD2G04.COM"
46
47 #include "CLM.h"
48
49 COMMON/SPEC2/KM,KINC,COEK,C3LAND(IO0,JM0),C3OICE(IO0,JM0) 5808.1
50 * ,C3LICE(IO0,JM0),WMGE(IO0,JM0),TSSFC(1,JM0,4) 5808.2
51 COMMON U,V,T,P,Q 5809.
52 COMMON/WORK1/CONV(IM0,JM0,LM0),PK(IM0,JM0,LM0),PREC(IM0,JM0),
53 & TPREC(IM0,JM0), 5810.
54 * COSZ1(IO0,JM0) 5811.
55 COMMON/WORK2/UT(IM0,JM0,LM0),VT(IM0,JM0,LM0),DU1(IO0,JM0),
56 & DV1(IO0,JM0), 5812.
57 * RA(8),ID(8),UMS(8) 5813.
58 COMMON/WORK3/E0(IO0,JM0,4),E1(IO0,JM0,4),EVAPOR(IO0,JM0,4), 5814.
59 * TGRND(IO0,JM0,4) 5814.1
60 COMMON/RDATA/ROUGHL(IO0,JM0) 5815.
61 LOGICAL POLE,PRNT,HPRNT
62 common/conprn/HPRNT
63 ! common/TSUR/TSURFC(JM0,0:13),TSURFT(JM0),TSURFD(JM0),DTEMSR(JM0)
64 #include "TSRF.COM"
65 common/SURRAD/TRSURF(JM0,4),SRSURF(JM0,4)
66 c REAL*8 B,TGV,TKV,TSV0,TSV1,TSV 5818.
67 integer IQ1,IQ2,IQ3
68 COMMON/CWMG/WMGEA(JM0),NWMGEA(JM0),CHAVER(JM0),DTAV(JM0),DQAV(JM0)
69 & ,Z0AV(JM0),WSAV(JM0),WS0AV(JM0),TAUAV(JM0)
70 C
71 COMMON/SURFLAND/ DUL1(JM0),DVL1(JM0),DT1L(JM0),DQ1L(JM0),
72 & WSSL(JM0),T2ML(JM0),
73 & TSSL(JM0),QSSL(JM0),USSL(JM0),VSSL(JM0),TAUSL(JM0),BLJ(JM0,50)
74 & ,ELHTG(JM0),SHTG(JM0),TAUXG(JM0),TAUYG(JM0)
75 c
76 DATA RVAP/461.5/ 5819.
77 DATA SHV/0./,SHW/4185./,SHI/2060./,RHOW/1000./,RHOI/916.6/, 5820.
78 * ALAMI/2.1762/,STBO/.5672573E-7/,TF/273.16/,TFO/-1.56/ 5821.
79 DATA Z1I/.1/,Z2LI/2.9/,Z1E/.1/,Z2E/4./,RHOS/91.66/,ALAMS/.35/ 5822.
80 QSAT(TM,PR,QLH)=3.797915*EXP(QLH*(7.93252E-6-2.166847E-3/TM))/PR 5836.
81 DLQSDT(TM,QLH)=QLH*2.166847E-3/(TM*TM)
82 DATA IFIRST/1/ 5838.
83 ROSNOW(X)=0.54*X/LOG(1.+0.54*X/275.)
84 ALSNOW(X)=2.8E-6*X**2
85 C**** 5839.
86 C**** FDATA 2 LAND COVERAGE (1) 5840.
87 C**** 3 RATIO OF LAND ICE COVERAGE TO LAND COVERAGE (1) 5841.
88 C**** 5842.
89 C**** ODATA 1 OCEAN TEMPERATURE (C) 5843.
90 C**** 2 RATIO OF OCEAN ICE COVERAGE TO WATER COVERAGE (1) 5844.
91 C**** 3 OCEAN ICE AMOUNT OF SECOND LAYER (KG/M**2) 5845.
92 C**** 5846.
93 C**** GDATA 1 OCEAN ICE SNOW AMOUNT (KG/M**2) 5847.
94 C**** 2 EARTH SNOW AMOUNT (KG/M**2) 5848.
95 C**** 3 OCEAN ICE TEMPERATURE OF FIRST LAYER (C) 5849.
96 C**** 4 EARTH TEMPERATURE OF FIRST LAYER (C) 5850.
97 C**** 5 EARTH WATER OF FIRST LAYER (KG/M**2) 5851.
98 C**** 6 EARTH ICE OF FIRST LAYER (KG/M**2) 5852.
99 C**** 7 OCEAN ICE TEMPERATURE OF SECOND LAYER (C) 5853.
100 C**** 8 EARTH TEMPERATURE OF SECOND LAYER (C) 5854.
101 C**** 9 EARTH WATER OF SECOND LAYER (KG/M**2) 5855.
102 C**** 10 EARTH ICE OF SECOND LAYER (KG/M**2) 5856.
103 C**** 12 LAND ICE SNOW AMOUNT (KG/M**2) 5857.
104 C**** 13 LAND ICE TEMPERATURE OF FIRST LAYER (C) 5858.
105 C**** 14 LAND ICE TEMPERATURE OF SECOND LAYER (C) 5859.
106 C**** 5860.
107 C**** BLDATA 1 COMPOSITE SURFACE WIND MAGNITUDE (M/S) 5861.
108 C**** 2 COMPOSITE SURFACE AIR TEMPERATURE (K) 5862.
109 C**** 3 COMPOSITE SURFACE AIR SPECIFIC HUMIDITY (1) 5863.
110 C**** 4 LAYER TO WHICH DRY CONVECTION MIXES (1) 5864.
111 C**** 5 FREE 5865.
112 C**** 6 COMPOSITE SURFACE U WIND 5866.
113 C**** 7 COMPOSITE SURFACE V WIND 5867.
114 C**** 8 COMPOSITE SURFACE MOMENTUM TRANSFER (TAU) 5868.
115 C**** 5869.
116 C**** VDATA 9 WATER FIELD CAPACITY OF FIRST LAYER (KG/M**2) 5870.
117 C**** 10 WATER FIELD CAPACITY OF SECOND LAYER (KG/M**2) 5871.
118 C**** 5872.
119 C**** 5874.
120 save
121 c print *,'sur4clm TAU=',TAU
122 NSTEPS=NSURF*NSTEP/NDYN 5875.
123 IF(IFIRST.NE.1) GO TO 50 5876.
124 print *,' SURFACE FOR CLM'
125 print *,' ZGS=30 m for LAND '
126 WMGMIN=8.
127 print *,'WMGMIN 4 LAND=',WMGMIN
128 IFIRST=0 5877.
129 print *,' WMGE'
130 print 258,(WMGE(1,J)+WMGMIN,J=1,JM)
131 258 format(12f5.1)
132 C SRCORX=1. 5878.13
133 CIAX=0.3
134 print *,' surfacen CIAX=',CIAX
135 print *,' QS=Q1, TS=T1'
136 print *,' WS=sqrt(0.75*W1+WGEM) '
137 IQ1=IM/4+1 5881.
138 IQ2=IM/2+1 5882.
139 IQ3=3*IM/4+1 5883.
140 ! DTSURF=NDYN*DT/NSURF 5884.
141 ! print *,' From SRF4CLM DTSURF=',DTSURF
142 ! DTSRCE=DT*NDYN 5885.
143 SHA=RGAS/KAPA 5886.
144 RVX=0. 5887.
145 50 CONTINUE
146 C**** ZERO OUT ENERGY AND EVAPORATION FOR GROUND AND INITIALIZE TGRND 5906.
147 DO 70 J=1,JM 5907.
148 DO 70 I=1,IM 5908.
149 TGRND(I,J,3)=GDATA(I,J,13) 5910.
150 TGRND(I,J,4)=GDATA(I,J,4) 5911.
151 70 CONTINUE
152 c print *,'After 70'
153 IHOUR=1.5+TOFDAY 5914.
154 C**** 5915.
155 C**** OUTSIDE LOOP OVER TIME STEPS, EXECUTED NSURF TIMES EVERY HOUR 5916.
156 C**** 5917.
157 C**** 5927.
158 C**** OUTSIDE LOOP OVER J AND I, EXECUTED ONCE FOR EACH GRID POINT 5928.
159 C**** 5929.
160 JPR=-7
161 DO 7000 J=1,JM 5930.
162 c print *,'After 7000 J=',J
163 c print *,IQ3
164 PRNT=.FALSE.
165 HEMI=1. 5931.
166 IF(J.LE.JM/2) HEMI=-1. 5932.
167 IF(J.EQ.1) GO TO 80 5936.
168 IF(J.EQ.JM) GO TO 90 5937.
169 WMG0=.5*(WMGE(1,J)+WMGE(1,J+1))+.001 5937.5
170 POLE=.FALSE. 5938.
171 IMAX=IM 5939.
172 GO TO 100 5940.
173 C**** CONDITIONS AT THE SOUTH POLE 5941.
174 80 POLE=.TRUE. 5942.
175 c print *,'After 80'
176 c print *,IQ1,IQ2,IQ3
177 IMAX=1 5943.
178 JVPO=2 5944.
179 RAPO=2.*RAPVN(1) 5945.
180 c print *,' RAPO=', RAPO
181 c II1=IQ1
182 c II2=IQ2
183 c II3=IQ3
184 c print *,II1,II2,II3
185 c print *,' III=',III
186 c print *,' U(IQ3,2,1)=',U(IQ3,2,1)
187 U1=.25*(U(1,2,1)+V(IQ1,2,1)-U(IQ2,2,1)-V(IQ3,2,1)) 5946.
188 V1=.25*(V(1,2,1)-U(IQ1,2,1)-V(IQ2,2,1)+U(IQ3,2,1)) 5947.
189 WMG0=WMGE(1,2) 5947.5
190 GO TO 100 5948.
191 C**** CONDITIONS AT THE NORTH POLE 5949.
192 90 POLE=.TRUE. 5950.
193 IMAX=1 5951.
194 JVPO=JM 5952.
195 RAPO=2.*RAPVS(JM) 5953.
196 U1=.25*(U(1,JM,1)-V(IQ1,JM,1)-U(IQ2,JM,1)+V(IQ3,JM,1)) 5954.
197 V1=.25*(V(1,JM,1)+U(IQ1,JM,1)-V(IQ2,JM,1)-U(IQ3,JM,1)) 5955.
198 WMG0=WMGE(1,JM) 5955.5
199 C**** ZERO OUT SURFACE DIAGNOSTICS WHICH WILL BE SUMMED OVER LONGITUDE 5956.
200 100 CONTINUE
201 c print *,'After 100'
202 BTS=0.
203 BWS=0.
204 BWMG=0.
205 IM1=IM 5969.
206 i=1
207 tsl4clm(i,j)=0.0
208 qs4clm(i,j)=0.0
209 ps4clm(i,j)=0.0
210 ws4clm(i,j)=0.0
211 us4clm(i,j)=0.0
212 vs4clm(i,j)=0.0
213 DO 6000 I=1,IMAX 5970.
214 C**** 5971.
215 C**** DETERMINE SURFACE CONDITIONS 5972.
216 C**** 5973.
217 PLAND=FDATA(I,J,2) 5974.
218 PWATER=1.-PLAND 5975.
219 PLICE=FDATA(I,J,3)*PLAND 5976.
220 PEARTH=PLAND-PLICE 5977.
221 POICE=ODATA(I,J,2)*PWATER 5978.
222 POCEAN=PWATER-POICE 5979.
223 if(POCEAN.LE.1.E-5)then
224 POCEAN=0.
225 POICE=PWATER
226 endif
227 TTOFR=PEARTH+PLICE+POICE+POCEAN
228 if(abs(TTOFR-1).gt.1.e-3)then
229 print *,' From surface TTOFR=',TTOFR
230 print *,' J=',J,' PLAND=',PLAND,' POCEAN=',POCEAN
231 print *,'POICE=',POICE,' ODATA(I,J,2)=',ODATA(I,J,2)
232 stop
233 end if
234 SP=P(I,J) 5980.
235 PS=SP+PTOP 5981.
236 PSK=EXPBYK(PS) 5982.
237 P1=SIG(1)*SP+PTOP 5983.
238 P1K=EXPBYK(P1) 5984.
239 IF(POLE) GO TO 1200 5993.
240 U1=.25*(U(IM1,J,1)+U(I,J,1)+U(IM1,J+1,1)+U(I,J+1,1)) 5994.
241 V1=.25*(V(IM1,J,1)+V(I,J,1)+V(IM1,J+1,1)+V(I,J+1,1)) 5995.
242 1200 TH1=T(I,J,1) 5996.
243 Q1=Q(I,J,1) 5997.
244 THV1=TH1*(1.+Q1*RVX) 5998.
245 C**** ZERO OUT QUANTITIES TO BE SUMMED OVER SURFACE TYPES 6002.
246 USS=0. 6003.
247 VSS=0. 6004.
248 WSS=0. 6005.
249 TSS=0. 6006.
250 QSS=0. 6007.
251 C**** 6032.
252 2400 IF(PLAND.LE.0.) GO TO 5000 6074.
253 ZGS=30. 6078.
254 IF(PLICE.LE.0.) GO TO 2600 6080.
255 C**** 6081.
256 C**** LAND ICE 6082.
257 C**** 6083.
258 ITYPE=3 6084.
259 PTYPE=PLICE 6085.
260 TG1=TGRND(I,J,3) 6087.
261 c ELHX=LHS 6094.
262 if (TG1.le.0.0)ELHX=LHS
263 if (TG1.gt.0.0)ELHX=LHE
264 GO TO 3000 6095.
265 C**** 6096.
266 2600 IF(PEARTH.LE.0.) GO TO 5000 6097.
267 C**** 6098.
268 C**** EARTH 6099.
269 C**** 6100.
270 ITYPE=4 6101.
271 PTYPE=PEARTH 6102.
272 TG1=TGRND(I,J,4) 6104.
273 if (TG1.le.0.0)ELHX=LHS
274 if (TG1.gt.0.0)ELHX=LHE
275 C**** 6134.
276 C**** BOUNDARY LAYER INTERACTION 6135.
277 C**** 6136.
278 3000 continue
279 if(J.eq.JPR)then
280 print *,' after 3000'
281 print *,'TAU=',TAU,' NS=',NS,' ITYPE=',ITYPE
282 print *,'CDH=',CDH,' RGAS=',RGAS
283 print *,'U1=',U1,' V1=',V1
284 print *,'WMGO=',WMGO
285 endif
286 TKV=THV1*PSK 6137.
287 C**** LOOP OVER GROUND TIME STEPS 6148.
288 TG=TG1+TF 6150.
289 QG=QSAT(TG,PS,ELHX) 6151.
290 TGV=TG*(1.+QG*RVX) 6152.
291 UG=0.75*U1
292 VG=0.75*V1
293 W1=SQRT(UG*UG+VG*VG)
294 WMG=WMG0+WMGMIN
295 WS=SQRT(W1**2+WMG)
296 RW=1.0
297 if(W1.ne.0.0)RW=WS/W1
298 THS=TH1
299 QS=Q1
300 TSV=THS*PSK
301 RIGS=ZGS*GRAV*(TSV-TGV)/(TGV*WS*WS)
302 IF(RIGS.LE.0) THEN
303 C surface layer has unstable stratification
304 CIA=TWOPI*0.0625/(1.+WS*CIAX)
305 ELSE
306 C surface layer has stable stratification
307 CIA=TWOPI*(0.09375-0.03125/(1.+4*RIGS**2))/(1.+WS*CIAX)
308 END IF
309 if(J.eq.JPR)then
310 print *,'TAU=',TAU,' NS=',NS,' ITYPE=',ITYPE
311 print *,'WS=',WS,' ZGS=',ZGS
312 print *,'RIGS=',RIGS,' TGV=',TGV
313 endif
314 USR=COS(CIA)
315 VSR=SIN(CIA)*HEMI
316 US=(USR*UG-VSR*VG)
317 VS=(VSR*UG+USR*VG)
318 if(J.eq.JPR)then
319 c print *,' '
320 print *,'TAU=',TAU,' NS=',NS,' ITYPE=',ITYPE
321 print *,'CDH=',CDH,' RGAS=',RGAS
322 print *,'PS=',PS,' TSV=',TSV
323 print *,'WS=',WS
324 endif
325 TS=TSV/(1.+QS*RVX) 6467.
326 c QSATS=QSAT(TS,PS,ELHX) 6468.
327 c IF(QS.LE.QSATS) GO TO 3500 6469.
328 c DQSSDT=QSATS*ELHX/(RVAP*TS*TS) 6470.
329 c X=(QS-QSATS)/(DQSSDT+(SHA/ELHX)) 6471.
330 c TS=TS+X 6472.
331 c QS=QS+X*(SHA/ELHX) 6473.
332 c3500 CONTINUE
333
334 if(ITYPE.EQ.4.or.ITYPE.EQ.3)then
335 tsl4clm(i,j)=tsl4clm(i,j)+TS*PTYPE/PLAND
336 qs4clm(i,j)=qs4clm(i,j)+QS*PTYPE/PLAND
337 ps4clm(i,j)=ps4clm(i,j)+PS*PTYPE/PLAND
338 ws4clm(i,j)=ws4clm(i,j)+WS*PTYPE/PLAND
339 ! us4clm(i,j)=us4clm(i,j)+US*PTYPE/PLAND
340 ! vs4clm(i,j)=vs4clm(i,j)+VS*PTYPE/PLAND
341 us4clm(i,j)=us4clm(i,j)+RW*US*PTYPE/PLAND
342 vs4clm(i,j)=vs4clm(i,j)+RW*VS*PTYPE/PLAND
343 endif
344
345 TSSFC(I,J,ITYPE)=TS 6521.5
346 USS=USS+US*PTYPE 6524.
347 VSS=VSS+VS*PTYPE 6525.
348 WSS=WSS+WS*PTYPE 6526.
349 TSS=TSS+TS*PTYPE 6527.
350 QSS=QSS+QS*PTYPE 6528.
351 GO TO (5000,5000,4400,4600),ITYPE 6551.
352 C**** 6552.
353 C**** LAND ICE 6569.
354 C**** 6570.
355 4400 CONTINUE
356 BTS=BTS+(TS-TF)*PLICE 6574.
357 BWS=BWS+WS*PLICE
358 BWMG=BWMG+SQRT(WMG)*PLICE
359 GO TO 2600 6575.
360 C**** 6576.
361 C**** EARTH 6577.
362 C**** 6578.
363 4600 CONTINUE
364 BTS=BTS+(TS-TF)*PEARTH 6582.
365 BWS=BWS+WS*PEARTH
366 BWMG=BWMG+SQRT(WMG)*PEARTH
367
368 C**** NON-OCEAN POINTS WHICH ARE NOT MELTING OR FREEZING WATER USE 6583.
369 C**** IMPLICIT TIME STEPS 6584.
370 C**** 6585.
371 C**** UPDATE SURFACE AND FIRST LAYER QUANTITIES 6586.
372 C**** 6587.
373 5000 CONTINUE
374 C**** 6596.
375 C**** ACCUMULATE DIAGNOSTICS 6597.
376 C**** 6598.
377 6000 IM1=I 6662.
378 WSSL(J)=WSS
379 TSSL(J)=TSS
380 QSSL(J)=QSS
381 USSL(J)=USS
382 VSSL(J)=VSS
383 C**** QUANTITIES ACCUMULATED FOR SURFACE TYPE TABLES IN DIAG1 6663.
384 BLJ(J,37)=BWS
385 BLJ(J,28)=BWMG
386 BLJ(J,23)=BTS
387 7000 CONTINUE 6677.
388 ! print *,' From SRF4CLM TAU=',TAU
389 ! i=1
390 ! j=1
391 ! print *,'TS(1),TS(2)'
392 ! print *,tsl4clm(1,1),tsl4clm(1,2)
393 ! print *,dsw4clm(i,j),dlw4clm(i,j)
394 ! print *,swinr4clm(i,j), swvis4clm(i,j)
395 C**** 6678.
396 #endif
397 RETURN 6795.
398 END 6824.

  ViewVC Help
Powered by ViewVC 1.1.22