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

Contents of /MITgcm_contrib/jscott/igsm/src/surf_land.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: +20 -15 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 SURF_LAND
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 CPL_CHEM )
38 !
39 #include "chem_para"
40 #include "chem_com"
41 !
42 #endif
43
44 #include "BD2G04.COM"
45
46 #if ( defined CLM )
47 #include "CLM.h"
48 #endif
49
50 COMMON/SPEC2/KM,KINC,COEK,C3LAND(IO0,JM0),C3OICE(IO0,JM0) 5808.1
51 * ,C3LICE(IO0,JM0),WMGE(IO0,JM0),TSSFC(1,JM0,4) 5808.2
52 COMMON U,V,T,P,Q 5809.
53 COMMON/WORK1/CONV(IM0,JM0,LM0),PK(IM0,JM0,LM0),PREC(IM0,JM0),
54 & TPREC(IM0,JM0), 5810.
55 * COSZ1(IO0,JM0) 5811.
56 COMMON/WORK2/UT(IM0,JM0,LM0),VT(IM0,JM0,LM0),DU1(IO0,JM0),
57 & DV1(IO0,JM0), 5812.
58 * RA(8),ID(8),UMS(8) 5813.
59 COMMON/WORK3/E0(IO0,JM0,4),E1(IO0,JM0,4),EVAPOR(IO0,JM0,4), 5814.
60 * TGRND(IO0,JM0,4) 5814.1
61 COMMON/RDATA/ROUGHL(IO0,JM0) 5815.
62 DIMENSION SINI(72),COSI(72) 5816.
63 LOGICAL POLE,PRNT,HPRNT
64 common/conprn/HPRNT
65 ! common/TSUR/TSURFC(JM0,0:13),TSURFT(JM0),TSURFD(JM0),DTEMSR(JM0)
66 #include "TSRF.COM"
67 common/SURRAD/TRSURF(JM0,4),SRSURF(JM0,4)
68 c REAL*8 B,TGV,TKV,TSV0,TSV1,TSV 5818.
69 COMMON/CWMG/WMGEA(JM0),NWMGEA(JM0),CHAVER(JM0),DTAV(JM0),DQAV(JM0)
70 & ,Z0AV(JM0),WSAV(JM0),WS0AV(JM0),TAUAV(JM0)
71 C
72 COMMON/SURFLAND/ DUL1(JM0),DVL1(JM0),DT1L(JM0),DQ1L(JM0),
73 & WSSL(JM0),T2ML(JM0),
74 & TSSL(JM0),QSSL(JM0),USSL(JM0),VSSL(JM0),TAUSL(JM0),BLJ(JM0,50)
75 & ,ELHTG(JM0),SHTG(JM0),TAUXG(JM0),TAUYG(JM0)
76 c
77 DATA RVAP/461.5/ 5819.
78 DATA SHV/0./,SHW/4185./,SHI/2060./,RHOW/1000./,RHOI/916.6/, 5820.
79 * ALAMI/2.1762/,STBO/.5672573E-7/,TF/273.16/,TFO/-1.56/ 5821.
80 DATA Z1I/.1/,Z2LI/2.9/,Z1E/.1/,Z2E/4./,RHOS/91.66/,ALAMS/.35/ 5822.
81 DIMENSION AROUGH(20),BROUGH(20),CROUGH(20),DROUGH(20),EROUGH(20) 5823.
82 DATA AROUGH/16.59,13.99,10.4,7.35,5.241,3.926,3.126,2.632,2.319, 5824.
83 *2.116,1.982,1.893,1.832,1.788,1.757,1.733,1.714,1.699,1.687,1.677/5825.
84 DATA BROUGH/3.245,1.733,0.8481,0.3899,0.1832,0.9026E-1,0.4622E-1, 5826.
85 * .241E-1,.1254E-1,.6414E-2,.3199E-2,.1549E-2,.7275E-3,.3319E-3, 5827.
86 * .1474E-3,.6392E-4,.2713E-4,.1130E-4,.4630E-5,.1868E-5/ 5828.
87 DATA CROUGH/5.111,3.088,1.682,.9239,.5626,.3994,.3282,.3017,.299 5829.
88 *,.3114,.3324,.3587,.3881,.4186,.4492,.4792,.5082,.5361,.5627, 5830.
89 * .5882/ 5831.
90 DATA DROUGH/1.24,1.02,0.806,0.682,0.661,0.771,0.797,0.895,0.994, 5832.
91 * 1.09,1.18,1.27,1.35,1.43,1.50,1.58,1.65,1.71,1.78,1.84/ 5833.
92 DATA EROUGH/0.128,0.130,0.141,0.174,0.238,0.330,0.438,0.550,0.660,5834.
93 * 0.766,0.866,0.962,1.05,1.14,1.22,1.30,1.37,1.45,1.52,1.58/ 5835.
94 QSAT(TM,PR,QLH)=3.797915*EXP(QLH*(7.93252E-6-2.166847E-3/TM))/PR 5836.
95 DLQSDT(TM,QLH)=QLH*2.166847E-3/(TM*TM)
96 c TLOG(Z0)=ALOG(.36*SQRTT/(FMAG*Z0))+2.302585*ROUGH-.08 5837.
97 DATA IFIRST/1/ 5838.
98 ROSNOW(X)=0.54*X/LOG(1.+0.54*X/275.)
99 ALSNOW(X)=2.8E-6*X**2
100 C**** 5839.
101 C**** FDATA 2 LAND COVERAGE (1) 5840.
102 C**** 3 RATIO OF LAND ICE COVERAGE TO LAND COVERAGE (1) 5841.
103 C**** 5842.
104 C**** ODATA 1 OCEAN TEMPERATURE (C) 5843.
105 C**** 2 RATIO OF OCEAN ICE COVERAGE TO WATER COVERAGE (1) 5844.
106 C**** 3 OCEAN ICE AMOUNT OF SECOND LAYER (KG/M**2) 5845.
107 C**** 5846.
108 C**** GDATA 1 OCEAN ICE SNOW AMOUNT (KG/M**2) 5847.
109 C**** 2 EARTH SNOW AMOUNT (KG/M**2) 5848.
110 C**** 3 OCEAN ICE TEMPERATURE OF FIRST LAYER (C) 5849.
111 C**** 4 EARTH TEMPERATURE OF FIRST LAYER (C) 5850.
112 C**** 5 EARTH WATER OF FIRST LAYER (KG/M**2) 5851.
113 C**** 6 EARTH ICE OF FIRST LAYER (KG/M**2) 5852.
114 C**** 7 OCEAN ICE TEMPERATURE OF SECOND LAYER (C) 5853.
115 C**** 8 EARTH TEMPERATURE OF SECOND LAYER (C) 5854.
116 C**** 9 EARTH WATER OF SECOND LAYER (KG/M**2) 5855.
117 C**** 10 EARTH ICE OF SECOND LAYER (KG/M**2) 5856.
118 C**** 12 LAND ICE SNOW AMOUNT (KG/M**2) 5857.
119 C**** 13 LAND ICE TEMPERATURE OF FIRST LAYER (C) 5858.
120 C**** 14 LAND ICE TEMPERATURE OF SECOND LAYER (C) 5859.
121 C**** 5860.
122 C**** BLDATA 1 COMPOSITE SURFACE WIND MAGNITUDE (M/S) 5861.
123 C**** 2 COMPOSITE SURFACE AIR TEMPERATURE (K) 5862.
124 C**** 3 COMPOSITE SURFACE AIR SPECIFIC HUMIDITY (1) 5863.
125 C**** 4 LAYER TO WHICH DRY CONVECTION MIXES (1) 5864.
126 C**** 5 FREE 5865.
127 C**** 6 COMPOSITE SURFACE U WIND 5866.
128 C**** 7 COMPOSITE SURFACE V WIND 5867.
129 C**** 8 COMPOSITE SURFACE MOMENTUM TRANSFER (TAU) 5868.
130 C**** 5869.
131 C**** VDATA 9 WATER FIELD CAPACITY OF FIRST LAYER (KG/M**2) 5870.
132 C**** 10 WATER FIELD CAPACITY OF SECOND LAYER (KG/M**2) 5871.
133 C**** 5872.
134 C**** ROUGHL LOG(ZGS/ROUGHNESS LENGTH) (LOGARITHM TO BASE 10) 5873.
135 C**** ROUGHL will be ROUGHNESS LENGTH
136 C**** 5874.
137 c print *,'surface TAU=',TAU
138 NSTEPS=NSURF*NSTEP/NDYN 5875.
139 IF(IFIRST.NE.1) GO TO 50 5876.
140 print *,' SURFACE FOR LAND ONLY'
141 print *,' ZGS=30 m for LAND '
142 WMGMIN=8.
143 WMGMIN=5.
144 print *,'WMGMIN 4 LAND=',WMGMIN
145 print *,'over land WMG=max(WMG0,WMGMIN)'
146 IFIRST=0 5877.
147 print *,' WMGE'
148 print 258,(WMGE(1,J),J=1,JM)
149 258 format(12f5.1)
150 ! print *,'ODATA(1,7,2)=',ODATA(1,7,2)
151 COEFSN=100./ROSNOW(10.)
152 COEFSN=1.
153 print *,' COEFSN=',COEFSN
154 do 2567 J=1,JM
155 NWMGEA(J)=0
156 WMGEA(J)=0.
157 CHAVER(J)=0.
158 DTAV(J)=0.
159 DQAV(J)=0.
160 Z0AV(J)=0.
161 WSAV(J)=0.
162 WS0AV(J)=0.
163 TAUAV(J)=0.
164 2567 CONTINUE
165 READ (519) ((ROUGHL(I,J),I=1,IO),J=1,JM) 5878.
166 DO 10 J=1,JM
167 ILAND=0.
168 SUM1=0. 5878.02
169 CONT1=0. 5878.03
170 CONT2=0.
171 DO 11 I=1,IO 5878.04
172 PLAND=C3LAND(I,J) 5878.05
173 CONT1=CONT1+PLAND 5878.06
174 ROUGHL(I,J)=10**(log10(30.)-ROUGHL(I,J))
175 C**** ROUGHL IS NOW ROUGHNESS LENGTH
176 11 SUM1=SUM1+PLAND*ROUGHL(I,J) 5878.07
177 IF(CONT1.LE.0.) GO TO 10 5878.08
178 SUM1=SUM1/CONT1 5878.09
179 DO 12 I=1,IO 5878.1
180 12 ROUGHL(I,J)=SUM1 5878.11
181 10 CONTINUE 5878.12
182 C SRCORX=1. 5878.13
183 CIAX=0.3
184 print *,' surfacen CIAX=',CIAX
185 print *,' QS=Q1, TS=T1'
186 print *,' WS=sqrt(0.75*W1+WGEM) '
187 print *,' ROUGHL'
188 print *,(ROUGHL(1,J),J=1,jm)
189 REWIND 519 5879.
190 LBLMM1=LBLM-1 5880.
191 IQ1=IM/4+1 5881.
192 IQ2=IM/2+1 5882.
193 IQ3=3*IM/4+1 5883.
194 DTSURF=NDYN*DT/NSURF 5884.
195 print *,' DTSURF=',DTSURF
196 DTSRCE=DT*NDYN 5885.
197 SHA=RGAS/KAPA 5886.
198 RVX=0. 5887.
199 ACE1I=Z1I*RHOI 5888.
200 HC1I=ACE1I*SHI 5889.
201 HC2LI=Z2LI*RHOI*SHI 5890.
202 HC1DE=Z1E*1129950. 5891.
203 HC2DE=Z2E*1129950.+3.5*.125*RHOW*3100. 5892.
204 Z1IBYL=Z1I/ALAMI 5893.
205 Z2LI3L=Z2LI/(3.*ALAMI) 5894.
206 BYRSL=1./(RHOS*ALAMS) 5895.
207 ZS1CO=.5*DSIG(1)*RGAS/GRAV 5896.
208 P1000K=EXPBYK(1000.) 5897.
209 COEFS=GRAV/(100.*DSIG(1)) 5898.
210 COEF1=(1.-SIG(2))/DSIGO(1) 5899.
211 COEF2=(SIG(1)-1.)/DSIGO(1) 5900.
212 DO 20 I=1,IM 5901.
213 SINI(I)=SIN((I-1)*TWOPI/FIM) 5902.
214 20 COSI(I)=COS((I-1)*TWOPI/FIM) 5903.
215 50 S0=S0X*1367./RSDIST 5904.
216 BYS0=RSDIST/1367. 5905.
217 C**** ZERO OUT ENERGY AND EVAPORATION FOR GROUND AND INITIALIZE TGRND 5906.
218 DO 70 J=1,JM 5907.
219 DO 70 I=1,IM 5908.
220 c TGRND(I,J,2)=GDATA(I,J,3) 5909.
221 TGRND(I,J,3)=GDATA(I,J,13) 5910.
222 TGRND(I,J,4)=GDATA(I,J,4) 5911.
223 DO 70 K=3,4 5912.
224 EVAPOR(I,J,K)=0.
225 E1(I,J,K)=0. 5913.
226 E0(I,J,K)=0. 5913.
227 70 CONTINUE
228 IHOUR=1.5+TOFDAY 5914.
229 C**** 5915.
230 C**** OUTSIDE LOOP OVER TIME STEPS, EXECUTED NSURF TIMES EVERY HOUR 5916.
231 C**** 5917.
232 DO 9000 NS=1,NSURF 5918.
233 MODDSF=MOD(NSTEPS+NS-1,NDASF) 5919.
234 C IF(MODDSF.EQ.0) IDACC(3)=IDACC(3)+1 5920.
235 MODD6=MOD(IDAY+NS,NSURF) 5921.
236 C**** ZERO OUT LAYER 1 WIND INCREMENTS 5922.
237 DO 60 J=1,JM 5923.
238 DUL1(J)=0.
239 60 DVL1(J)=0.
240 C**** 5927.
241 C**** OUTSIDE LOOP OVER J AND I, EXECUTED ONCE FOR EACH GRID POINT 5928.
242 C**** 5929.
243 JPR=-7
244 DO 7000 J=1,JM 5930.
245 ELHTG(J)=0.0
246 SHTG(J)=0.0
247 TAUXG(J)=0.0
248 TAUYG(J)=0.0
249 PRNT=j.eq.8
250 PRNT=.FALSE.
251 if(PRNT)then
252 if(ns.eq.1)then
253 write(78,*) ,' '
254 write(78,*) ,'TAU=',TAU
255 endif
256 write(78,*),'NS=',ns
257 endif
258 HEMI=1. 5931.
259 IF(J.LE.JM/2) HEMI=-1. 5932.
260 FCOR=2.*OMEGA*SINP(J) 5933.
261 FMAG=FCOR*HEMI 5934.
262 ROOT2F=SQRT(2.*FMAG) 5935.
263 IF(J.EQ.1) GO TO 80 5936.
264 IF(J.EQ.JM) GO TO 90 5937.
265 WMG0=.5*(WMGE(1,J)+WMGE(1,J+1))+.001 5937.5
266 POLE=.FALSE. 5938.
267 IMAX=IM 5939.
268 GO TO 100 5940.
269 C**** CONDITIONS AT THE SOUTH POLE 5941.
270 80 POLE=.TRUE. 5942.
271 IMAX=1 5943.
272 JVPO=2 5944.
273 RAPO=2.*RAPVN(1) 5945.
274 U1=.25*(U(1,2,1)+V(IQ1,2,1)-U(IQ2,2,1)-V(IQ3,2,1)) 5946.
275 V1=.25*(V(1,2,1)-U(IQ1,2,1)-V(IQ2,2,1)+U(IQ3,2,1)) 5947.
276 WMG0=WMGE(1,2)+.001 5947.5
277 GO TO 100 5948.
278 C**** CONDITIONS AT THE NORTH POLE 5949.
279 90 POLE=.TRUE. 5950.
280 IMAX=1 5951.
281 JVPO=JM 5952.
282 RAPO=2.*RAPVS(JM) 5953.
283 U1=.25*(U(1,JM,1)-V(IQ1,JM,1)-U(IQ2,JM,1)+V(IQ3,JM,1)) 5954.
284 V1=.25*(V(1,JM,1)+U(IQ1,JM,1)-V(IQ2,JM,1)-U(IQ3,JM,1)) 5955.
285 WMG0=WMGE(1,JM)+.001 5955.5
286 C**** ZERO OUT SURFACE DIAGNOSTICS WHICH WILL BE SUMMED OVER LONGITUDE 5956.
287 100 ATRHDT=0. 5957.
288 BTRHDT=0. 5958.
289 CTRHDT=0. 5959.
290 ASHDT=0. 5960.
291 BSHDT=0. 5961.
292 CSHDT=0. 5962.
293 AEVHDT=0. 5963.
294 BEVHDT=0. 5964.
295 CEVHDT=0. 5965.
296 ATS=0. 5966.
297 BTS=0. 5967.
298 CTS=0. 5968.
299 AT2=0. 5966.
300 BT2=0. 5967.
301 CT2=0. 5968.
302 ATAUL=0.
303 ATAUF=0.
304 BTAUL=0.
305 BTAUF=0.
306 CTAUL=0.
307 CTAUF=0.
308 AWS=0.
309 BWS=0.
310 CWS=0.
311 AWMG=0.
312 BWMG=0.
313 CWMG=0.
314 ACH=0.
315 BCH=0.
316 CCH=0.
317 IM1=IM 5969.
318 #if ( defined CLM )
319 if(NS.eq.1)then
320 i=1
321 tsl4clm(i,j)=0.0
322 qs4clm(i,j)=0.0
323 ps4clm(i,j)=0.0
324 ws4clm(i,j)=0.0
325 us4clm(i,j)=0.0
326 vs4clm(i,j)=0.0
327 endif
328 #endif
329 DO 6000 I=1,IMAX 5970.
330 C**** 5971.
331 C**** DETERMINE SURFACE CONDITIONS 5972.
332 C**** 5973.
333 PLAND=FDATA(I,J,2) 5974.
334 PWATER=1.-PLAND 5975.
335 PLICE=FDATA(I,J,3)*PLAND 5976.
336 PEARTH=PLAND-PLICE 5977.
337 POICE=ODATA(I,J,2)*PWATER 5978.
338 POCEAN=PWATER-POICE 5979.
339 if(POCEAN.LE.1.E-5)then
340 POCEAN=0.
341 POICE=PWATER
342 endif
343 TTOFR=PEARTH+PLICE+POICE+POCEAN
344 if(abs(TTOFR-1).gt.1.e-3)then
345 print *,' From surface TTOFR=',TTOFR
346 print *,' J=',J,' PLAND=',PLAND,' POCEAN=',POCEAN
347 print *,'POICE=',POICE,' ODATA(I,J,2)=',ODATA(I,J,2)
348 stop
349 end if
350 SP=P(I,J) 5980.
351 PS=SP+PTOP 5981.
352 PSK=EXPBYK(PS) 5982.
353 P1=SIG(1)*SP+PTOP 5983.
354 P1K=EXPBYK(P1) 5984.
355 WSOLD=BLDATA(I,J,1) 5985.
356 USOLD=BLDATA(I,J,6) 5986.
357 VSOLD=BLDATA(I,J,7) 5987.
358 TOLD=BLDATA(I,J,8) 5988.
359 SQRTT=SQRT(TOLD) 5989.
360 GKBYFW=.1296*GRAV/(FCOR*FMAG*WSOLD+1.E-20) 5990.
361 COSWS=GKBYFW*USOLD 5991.
362 SINWS=GKBYFW*VSOLD 5992.
363 IF(POLE) GO TO 1200 5993.
364 U1=.25*(U(IM1,J,1)+U(I,J,1)+U(IM1,J+1,1)+U(I,J+1,1)) 5994.
365 V1=.25*(V(IM1,J,1)+V(I,J,1)+V(IM1,J+1,1)+V(I,J+1,1)) 5995.
366 if(J.eq.JPR.or.J.eq.-12)then
367 print *,' J=',J
368 print *,'TAU=',TAU,' NS=',NS,' ITYPE=',ITYPE
369 print *,'U(I,J,1)=',U(I,J,1),' V(I,J,1)=',V(I,J,1)
370 print *,'U(I,J+1,1)=',U(I,J+1,1),' V(I,J+1,1)=',V(I,J+1,1)
371 print *,'U(IM1,J,1)=',U(IM1,J,1),' V(IM1,J,1)=',V(IM1,J,1)
372 print *,'U(IM1,J+1,1)=',U(IM1,J+1,1),
373 & ' V(IM1,J+1,1)=',V(IM1,J+1,1)
374 endif
375 1200 TH1=T(I,J,1) 5996.
376 Q1=Q(I,J,1) 5997.
377 DTH1=0.0
378 DQQ1=0.0
379 THV1=TH1*(1.+Q1*RVX) 5998.
380 c SRHEAT=SRHR(I,J,1)*COSZ1(I,J)*SRCOR 5999.
381 c SRHDT=SRHEAT*DTSURF 6000.
382 RMBYA=100.*SP*DSIG(1)/GRAV 6001.
383 C**** ZERO OUT QUANTITIES TO BE SUMMED OVER SURFACE TYPES 6002.
384 USS=0. 6003.
385 VSS=0. 6004.
386 WSS=0. 6005.
387 TSS=0. 6006.
388 T2MS=0.
389 QSS=0. 6007.
390 TAUS=0. 6008.
391 SINAPS=0. 6009.
392 COSAPS=0. 6010.
393 JR=J
394 DXYPJ=DXYP(J) 6012.
395 TG1S=0. 6013.
396 QGS=0. 6014.
397 BETAS=0. 6015.
398 TRHDTS=0. 6016.
399 SHDTS=0. 6017.
400 EVHDTS=0. 6018.
401 UGS=0. 6019.
402 VGS=0. 6020.
403 WGS=0. 6021.
404 USRS=0. 6022.
405 VSRS=0. 6023.
406 RIS1S=0. 6024.
407 RIGSS=0. 6025.
408 CDMS=0. 6026.
409 CDHS=0. 6027.
410 DGSS=0. 6028.
411 EDS1S=0. 6029.
412 PPBLS=0. 6030.
413 EVAPS=0. 6031.
414 C**** 6032.
415 2400 IF(PLAND.LE.0.) GO TO 5000 6074.
416 NGRNDZ=NGRND 6075.
417 ROUGH=ROUGHL(I,J) 6076.
418 ZGS=30. 6078.
419 ! WMGMIN=5.
420 TRHT0=TRSURF(J,2)
421 SRHEAT=SRSURF(J,2)*COSZ1(I,J)*SRCOR
422 IF(PLICE.LE.0.) GO TO 2600 6080.
423 C**** 6081.
424 C**** LAND ICE 6082.
425 C**** 6083.
426 ITYPE=3 6084.
427 PTYPE=PLICE 6085.
428 SNOW=GDATA(I,J,12) 6086.
429 TG1=TGRND(I,J,3) 6087.
430 TG2=GDATA(I,J,14) 6088.
431 if (SNOW.gt.10.)then
432 RHOS0=ROSNOW(SNOW)
433 else
434 RHOS0=275.
435 endif
436 RHOS=COEFSN*RHOS0
437 ALAMS=ALSNOW(RHOS0)
438 BYRSL=1./(RHOS*ALAMS)
439 c Z1BY6L=(Z1IBYL+SNOW*BYRSL)*.1666667 6089.
440 CDTERM=TG2 6090.
441 c CDENOM=1./(2.*Z1BY6L+Z2LI3L) 6091.
442 Z1BY2L=(Z1IBYL+SNOW*BYRSL)*0.5
443 CDENOM=1./(Z1BY2L+3.*Z2LI3L/2.)
444 HC1=HC1I+SNOW*SHI 6092.
445 BETA=1. 6093.
446 ELHX=LHS 6094.
447 GO TO 3000 6095.
448 C**** 6096.
449 2600 IF(PEARTH.LE.0.) GO TO 5000 6097.
450 C**** 6098.
451 C**** EARTH 6099.
452 C**** 6100.
453 ITYPE=4 6101.
454 PTYPE=PEARTH 6102.
455 SNOW=GDATA(I,J,2) 6103.
456 TG1=TGRND(I,J,4) 6104.
457 WTR1=GDATA(I,J,5) 6105.
458 ACE1=GDATA(I,J,6) 6106.
459 TG2=GDATA(I,J,8) 6107.
460 WTR2=GDATA(I,J,9) 6108.
461 ACE2=GDATA(I,J,10) 6109.
462 WFC1=VDATA(I,J,9) 6110.
463 WFC2=VDATA(I,J,10) 6111.
464 WTR1DRY=0.025*WFC1
465 HC1=HC1DE+WTR1*SHW+(ACE1+SNOW)*SHI 6112.
466 ALAM1D=2.+.5*(1.+2.*WTR1/WFC1) 6113.
467 ALAM2D=4. 6114.
468 RMULCH=1. 6115.
469 IF((SINP(J).GT..5).AND.(JDAY-91)*(273-JDAY).LT.0) RMULCH=.25 6116.
470 IF((SINP(J).LT.-.5).AND.(JDAY-91)*(273-JDAY).GE.0) RMULCH=.25 6117.
471 ALAM1V=RMULCH*(.4185+1.2555*WTR1/WFC1+ALAMI*ACE1/(Z1E*RHOI)) 6118.
472 ALAM3V=.8370 6119.
473 IF(TG2.LT.0.) ALAM3V=.4185+ALAMI*.15 6120.
474 ALAM2V=.125*(.4185+1.2555*WTR2/WFC2+ALAMI*ACE2/(5.*Z1E*RHOI)) 6121.
475 * +.875*ALAM3V 6122.
476 ALAM1E=VDATA(I,J,1)*ALAM1D+(1.-VDATA(I,J,1))*ALAM1V 6123.
477 ALAM2E=VDATA(I,J,1)*ALAM2D+(1.-VDATA(I,J,1))*ALAM2V 6124.
478 if (SNOW.gt.10.)then
479 RHOS0=ROSNOW(SNOW)
480 else
481 RHOS0=275.
482 endif
483 RHOS=COEFSN*RHOS0
484 ALAMS=ALSNOW(RHOS0)
485 BYRSL=1./(RHOS*ALAMS)
486 c Z1BY6L=(Z1E/ALAM1E+SNOW*BYRSL)*.1666667 6125.
487 Z1BY2L=(Z1E/ALAM1E+SNOW*BYRSL)*0.5
488 CDTERM=TG2 6126.
489 c CDENOM=1./(2.*Z1BY6L+Z2E/(3.*ALAM2E)) 6127.
490 CDENOM=1./(Z1BY2L+Z2E/(2.*ALAM2E))
491 BETA=1. 6128.
492 ELHX=LHS 6129.
493 IF(SNOW.GT.0.) GO TO 3000 6130.
494 BETA=(WTR1+ACE1)/WFC1 6131.
495 BETA=max(((WTR1+ACE1-WTR1DRY)/WFC1),0.0)
496 PFROZN=ACE1/(WTR1+ACE1+1.E-20) 6132.
497 ELHX=LHE+LHM*PFROZN 6133.
498 HC2E=HC2DE+WTR2*SHW+ACE2*SHI
499 C**** 6134.
500 C**** BOUNDARY LAYER INTERACTION 6135.
501 C**** 6136.
502 3000 continue
503 SRHDT=SRHEAT*DTSURF
504 TKV=THV1*PSK 6137.
505 ZS1=ZS1CO*TKV*SP/PS 6138.
506 P1=SIG(1)*SP+PTOP 6139.
507 DTGRND=DTSURF/NGRNDZ 6143.
508 SHDT=0. 6144.
509 EVHDT=0. 6145.
510 TRHDT=0. 6146.
511 F1DT=0. 6147.
512 C**** LOOP OVER GROUND TIME STEPS 6148.
513 DO 3600 NG=1,NGRNDZ 6149.
514 TG=TG1+TF 6150.
515 QG=QSAT(TG,PS,ELHX) 6151.
516 TGV=TG*(1.+QG*RVX) 6152.
517 ! UG=U1
518 ! VG=V1
519 UG=0.75*U1 ! 07/14/2006
520 VG=0.75*V1 ! 07/14/2006
521 W1=SQRT(U1*U1+V1*V1)
522 W1=SQRT(UG*UG+VG*VG) ! 07/14/2006
523 WS0=W1
524 ! WMG=WMG0+WMGMIN
525 ! 07/17/2006
526 WMG=WMG0+WMGMIN
527 WS=SQRT(W1**2+WMG)
528 if(J.eq.JPR)then
529 print *,' '
530 print *,'TAU=',TAU,' NS=',NS,' ITYPE=',ITYPE
531 print *,'TG=',TG,' QG=',QG
532 print *,'RVX=',RVX,' TG1=',TG1
533 endif
534
535 WG=WS
536 THS=TH1
537 QS=Q1
538 TSV=THS*PSK
539 Z0=ROUGH
540 ROUGH=log10(ZGS/ROUGH)
541 CDN=.0231/(ROUGH*ROUGH)
542 LR=ROUGH*2.-.5
543 IF(LR.GT.20) LR=20
544 IF(LR.LT.1) LR=1
545 RIGS=ZGS*GRAV*(TSV-TGV)/(TGV*WS*WS)
546 SINAP=0.
547 COSAP=1.
548 IF(RIGS.LE.0) THEN
549 C surface layer has unstable stratification
550 CIA=TWOPI*0.0625/(1.+WS*CIAX)
551 DM=SQRT((1.-AROUGH(LR)*RIGS)*(1.-BROUGH(LR)*RIGS)/
552 * (1.-CROUGH(LR)*RIGS))
553 DH=1.35*SQRT((1.-DROUGH(LR)*RIGS)/(1.-EROUGH(LR)*RIGS))
554 ELSE
555 C surface layer has stable stratification
556 CIA=TWOPI*(0.09375-0.03125/(1.+4*RIGS**2))/(1.+WS*CIAX)
557 DM=1./(1.+(11.238+89.9*RIGS)*RIGS)
558 DH=1.35/(1.+1.93*RIGS)
559 END IF
560 CDH=CDN*DM*DH
561 if(J.eq.JPR)then
562 print *,'TAU=',TAU,' NS=',NS,' ITYPE=',ITYPE
563 print *,'WS=',WS,' ZGS=',ZGS
564 print *,'DM=',DM,' DH=',DH
565 print *,'RIGS=',RIGS,' TGV=',TGV
566 endif
567 USR=COS(CIA)
568 VSR=SIN(CIA)*HEMI
569 US=(USR*UG-VSR*VG)
570 VS=(VSR*UG+USR*VG)
571 RCDHWS=CDH*WS*100.*PS/(RGAS*TSV)
572 if(J.eq.JPR)then
573 c print *,' '
574 print *,'TAU=',TAU,' NS=',NS,' ITYPE=',ITYPE
575 print *,'CDH=',CDH,' RGAS=',RGAS
576 print *,'PS=',PS,' TSV=',TSV
577 print *,'WS=',WS,' RCDHWS=',RCDHWS
578 endif
579 TS=TSV/(1.+QS*RVX) 6467.
580 QSATS=QSAT(TS,PS,ELHX) 6468.
581 IF(QS.LE.QSATS) GO TO 3500 6469.
582 DQSSDT=QSATS*ELHX/(RVAP*TS*TS) 6470.
583 X=(QS-QSATS)/(DQSSDT+(SHA/ELHX)) 6471.
584 TS=TS+X 6472.
585 QS=QS+X*(SHA/ELHX) 6473.
586 3500 CONTINUE
587
588 #if ( defined CLM )
589 if(ITYPE.EQ.4.or.ITYPE.EQ.3)then
590 tsl4clm(i,j)=tsl4clm(i,j)+TS*PTYPE/PLAND
591 qs4clm(i,j)=qs4clm(i,j)+QS*PTYPE/PLAND
592 ps4clm(i,j)=ps4clm(i,j)+PS*PTYPE/PLAND
593 ws4clm(i,j)=ws4clm(i,j)+WS*PTYPE/PLAND
594 us4clm(i,j)=us4clm(i,j)+US*PTYPE/PLAND
595 vs4clm(i,j)=vs4clm(i,j)+VS*PTYPE/PLAND
596 endif
597 #endif
598
599 C**** CALCULATE RHOS*CDM*WS AND RHOS*CDH*WS 6474.
600 CDM=CDN*DM 6475.
601 RCDMWS=CDM*WS*100.*PS/(RGAS*TS) 6476.
602 C**** CALCULATE FLUXES OF SENSIBLE HEAT, LATENT HEAT, THERMAL 6478.
603 C**** RADIATION, AND CONDUCTION HEAT (WATTS/M**2) 6479.
604 SHEAT=SHA*RCDHWS*(TS-TG) 6480.
605 BETAUP=BETA 6481.
606 IF(QS.GT.QG) BETAUP=1. 6482.
607 EVHEAT=(LHE+TG1*SHV)*BETAUP*RCDHWS*(QS-QG) 6483.
608 c TRHEAT=TRHR(I,J,1)-STBO*(TG*TG)*(TG*TG) 6484.
609 TRHEAT=TRHT0-STBO*(TG*TG)*(TG*TG)
610 if(J.eq.JPR)then
611 c print *,' '
612 print *,'TAU=',TAU,' NS=',NS,' ITYPE=',ITYPE
613 print *,'TRHT0=',TRHT0,' STBO=',STBO
614 print *,'TG=',TG,' TS=',TS
615 print *,'TRHEAT=',TRHEAT
616 print *,'SHA=',SHA,' RCDHWS=',RCDHWS
617 print *,'SHEAT=',SHEAT
618 endif
619 TG1OLD=TG1
620 SHEATOLD=SHEAT
621 C**** CALCULATE FLUXES USING IMPLICIT TIME STEP FOR NON-OCEAN POINTS 6486.
622 F0=SRHEAT+TRHEAT+SHEAT+EVHEAT 6487.
623 c F1=(TG1-CDTERM-F0*Z1BY6L)*CDENOM 6488.
624 F1=(TG1-CDTERM)*CDENOM
625 DSHDTG=-RCDHWS*SHA
626 DQGDTG=QG*ELHX/(RVAP*TG*TG) 6490.
627 DEVDTG=-RCDHWS*LHE*BETAUP*DQGDTG
628 DTRDTG=-4.*STBO*TG*TG*TG 6492.
629 DF0DTG=DSHDTG+DEVDTG+DTRDTG 6493.
630 c DFDTG=DF0DTG-(1.-DF0DTG*Z1BY6L)*CDENOM 6493.5
631 DFDTG=DF0DTG-CDENOM
632 c DF1DTG=(1.-DF0DTG*Z1BY6L)*CDENOM
633 DF1DTG=CDENOM
634 DTG=(F0-F1)*DTGRND/(HC1-DTGRND*DFDTG) 6494.
635 SHDT=SHDT+DTGRND*(SHEAT+DTG*DSHDTG) 6495.
636 EVHDT=EVHDT+DTGRND*(EVHEAT+DTG*DEVDTG) 6496.
637 TRHDT=TRHDT+DTGRND*(TRHEAT+DTG*DTRDTG) 6497.
638 TG1=TG1+DTG
639 c F1DT=F1DT+DTGRND*(TG1-CDTERM-(F0+DTG*DF0DTG)*Z1BY6L)*CDENOM 6498.
640 F1DT=F1DT+DTGRND*(TG1-CDTERM)*CDENOM
641 DUL1(J)=DUL1(J)+PTYPE*DTGRND*RCDMWS*US*COEFS/SP
642 DVL1(J)=DVL1(J)+PTYPE*DTGRND*RCDMWS*VS*COEFS/SP
643 c TG1=TG1+DTG 6501.
644 3600 CONTINUE 6502.
645 EPS=1.D-8
646 c print *,'FROM SURFACE NS=',NS
647 c print *,'J=',J,' ITYPE=',ITYPE
648 c print *,RCDMWS,WS
649 WWS=max(W1,1.D-4)
650 c RO=SP*100/(RGAS*TG)
651 c print *,'RO=',RO
652 c USTAR=SQRT(RCDMWS*WS/RO)
653 c TSTAR=SHEATOLD/(0.35*1007.*RO*USTAR)
654 c ALPHAH=DH
655 c TT2M=TG+TSTAR/ALPHAH*LOG(2.0/Z0)
656 c TT2M=TG+TSTAR/ALPHAH*LOG(ZGS/Z0)
657 c print *,'RIGS=',RIGS,' Z0=',Z0
658 c print *,'CDN=',CDN
659 c print *,'H=',SHDT/DTSURF,' TGM=',RCDMWS*WS
660 c print *,' SHEATOLD=',SHEATOLD
661 c print *,' USTAR=',USTAR,' TSTAR=',TSTAR
662 c print *,' ALPHAH=',ALPHAH,' TT2M=',TT2M
663 c print *,' TT2M=',TT2M
664 ZTEM=ZGS
665 ZTEM=2.0
666 c print *,'ZTEM=',ZTEM
667 CALL TZM(T2M,TH2M,ZTEM,Z0,ZGS,SP,TG,TS,RIGS,WS,
668 & -SHEATOLD,RCDMWS*WS,LR,EPS)
669 c print *,'FROM SURFACE'
670 c print *,'TS=',TS,' TG=',TG
671 c print *,' T2M=',T2M,' TH2M=',TH2M
672 F0DT=CORSR*SRHDT+TRHDT+SHDT+EVHDT 6510.
673 if(J.eq.JPR)then
674 print *,'TAU=',TAU,' NS=',NS,' ITYPE=',ITYPE
675 print *,'DTSURF=',DTSURF,' CORSR=',CORSR
676 print *,'SRHDT=',SRHDT,' TRHDT=',TRHDT
677 print *,'SHDT=',SHDT,' EVHDT=',EVHDT
678 print *,'F0DT=',F0DT
679 print *,'US=',US,' VS=',VS
680 print *,'COEFS=',COEFS,' SP=',SP
681 endif
682 c print *,'From surface ',TAU,CORSR,SRHDT,TRHDT,SHDT,EVHDT
683 C**** CALCULATE EVAPORATION 6511.
684 CCC DQ1=EVHDT/((LHE+TG1*SHV)*RMBYA) 6512.
685 DQ1=EVHDT/(ELHX*RMBYA)
686 IF(DQ1*PTYPE.LE.Q1) GO TO 3720 6513.
687 DQ1=Q1/PTYPE 6514.
688 CCC EVHDT=DQ1*(LHE+TG1*SHV)*RMBYA 6515.
689 EVHDT=DQ1*ELHX*RMBYA
690 3720 EVAP=-DQ1*RMBYA 6516.
691 C**** ACCUMULATE SURFACE FLUXES AND PROGNOSTIC AND DIAGNOSTIC QUANTITIES6517.
692 E0(I,J,ITYPE)=E0(I,J,ITYPE)+F0DT 6518.
693 E1(I,J,ITYPE)=E1(I,J,ITYPE)+F1DT 6519.
694 EVAPOR(I,J,ITYPE)=EVAPOR(I,J,ITYPE)+EVAP 6520.
695 if(PRNT)then
696 c write(78,*) ,' '
697 c write(78,*) ,'TAU=',TAU
698 write(78,*) ,'J=',j,' ITYPE=',ITYPE,' PTYPE=',PTYPE
699 write(78,*) ,'TS=',TS,' TG=',TG,' QS=',QS
700 write(78,*) ,'TG1=',TG1,' TG1OLD=',TG1OLD
701 write(78,*) ,'TG2=',TG2
702 write(78,*) ,'SHEAT=',SHEAT,' EVHEAT=',EVHEAT
703 write(78,*) ,'TRHEAT=',TRHEAT,' SRHEAT=',SRHEAT
704 write(78,*) ,'EVAP mm/day=',24.*3600.*EVAP/DTSURF
705 write(78,*) ,'EVAP=',EVAP,
706 & ' F0DT=',F0DT/DTSURF,' F1DT=',F1DT/DTSURF
707 endif
708 TGRND(I,J,ITYPE)=TG1 6521.
709 TSSFC(I,J,ITYPE)=TS 6521.5
710 c TH1=TH1-SHDT*PTYPE/(SHA*RMBYA*P1K) 6522.
711 c Q1=Q1-DQ1*PTYPE 6523.
712
713 DTH1=DTH1-SHDT*PTYPE/(SHA*RMBYA*P1K)
714 DQQ1=DQQ1-DQ1*PTYPE
715
716
717 USS=USS+US*PTYPE 6524.
718 VSS=VSS+VS*PTYPE 6525.
719 WSS=WSS+WS*PTYPE 6526.
720 TSS=TSS+TS*PTYPE 6527.
721 T2MS=T2MS+T2M*PTYPE
722 QSS=QSS+QS*PTYPE 6528.
723 ELHTG(J)=ELHTG(J)+EVHEAT*PTYPE
724 SHTG(J)=SHTG(J)+SHEAT*PTYPE
725 TAUXG(J)=TAUXG(J)+RCDMWS*US*PTYPE
726 TAUYG(J)=TAUYG(J)+RCDMWS*VS*PTYPE
727 ! TAUS=TAUS+CDM*WS*W1*PTYPE 6529.
728 TAUS=TAUS+CDM*WS*WS*PTYPE ! 7/14/2006
729 SINAPS=SINAPS+SINAP*PTYPE 6530.
730 COSAPS=COSAPS+COSAP*PTYPE 6531.
731 TG1S=TG1S+TG1*PTYPE 6532.
732 QGS=QGS+QG*PTYPE 6533.
733 BETAS=BETAS+BETA*PTYPE 6534.
734 TRHDTS=TRHDTS+TRHDT*PTYPE 6535.
735 SHDTS=SHDTS+SHDT*PTYPE 6536.
736 EVHDTS=EVHDTS+EVHDT*PTYPE 6537.
737 UGS=UGS+UG*PTYPE 6538.
738 VGS=VGS+VG*PTYPE 6539.
739 WGS=WGS+WG*PTYPE 6540.
740 USRS=USRS+USR*PTYPE 6541.
741 VSRS=VSRS+VSR*PTYPE 6542.
742 c RIS1S=RIS1S+RIS1*PTYPE 6543.
743 RIGSS=RIGSS+RIGS*PTYPE 6544.
744 CDMS=CDMS+CDM*PTYPE 6545.
745 CDHS=CDHS+CDH*PTYPE 6546.
746 c DGSS=DGSS+DGS*PTYPE 6547.
747 c EDS1S=EDS1S+EDS1*PTYPE 6548.
748 c PPBLS=PPBLS+PPBL*PTYPE 6549.
749 EVAPS=EVAPS+EVAP*PTYPE 6550.
750 GO TO (5000,5000,4400,4600),ITYPE 6551.
751 C**** 6552.
752 C**** LAND ICE 6569.
753 C**** 6570.
754 4400 BSHDT=BSHDT+SHDT*PLICE 6571.
755 BEVHDT=BEVHDT+EVHDT*PLICE 6572.
756 BTRHDT=BTRHDT+TRHDT*PLICE 6573.
757 BTS=BTS+(TS-TF)*PLICE 6574.
758 c BT2=BT2+(TH2M-TF)*PLICE
759 BT2=BT2+(T2M-TF)*PLICE
760 BTAUL=BTAUL+RCDMWS*US*PLICE
761 BTAUF=BTAUF+RCDMWS*VS*PLICE
762 BWS=BWS+WS*PLICE
763 BWMG=BWMG+SQRT(WMG)*PLICE
764 BCH=BCH+RCDHWS*PLICE
765 GO TO 2600 6575.
766 C**** 6576.
767 C**** EARTH 6577.
768 C**** 6578.
769 4600 BSHDT=BSHDT+SHDT*PEARTH 6579.
770 BEVHDT=BEVHDT+EVHDT*PEARTH 6580.
771 BTRHDT=BTRHDT+TRHDT*PEARTH 6581.
772 BTS=BTS+(TS-TF)*PEARTH 6582.
773 c BT2=BT2+(TH2M-TF)*PEARTH
774 BT2=BT2+(T2M-TF)*PEARTH
775 BTAUL=BTAUL+RCDMWS*US*PEARTH
776 BTAUF=BTAUF+RCDMWS*VS*PEARTH
777 BWS=BWS+WS*PEARTH
778 BWMG=BWMG+SQRT(WMG)*PEARTH
779 BCH=BCH+RCDHWS*PEARTH
780
781 C**** NON-OCEAN POINTS WHICH ARE NOT MELTING OR FREEZING WATER USE 6583.
782 C**** IMPLICIT TIME STEPS 6584.
783 C**** 6585.
784 C**** UPDATE SURFACE AND FIRST LAYER QUANTITIES 6586.
785 C**** 6587.
786 5000 CONTINUE
787 DT1L(J)=DTH1
788 DQ1L(J)=DQQ1
789 WSSL(J)=WSS
790 TSSL(J)=TSS
791 T2ML(J)=T2MS
792 QSSL(J)=QSS
793 USSL(J)=USS
794 VSSL(J)=VSS
795 TAUSL(J)=TAUS
796 ELHTG(J)=ELHTG(J)/PLAND
797 SHTG(J)=SHTG(J)/PLAND
798 TAUXG(J)=TAUXG(J)/PLAND
799 TAUYG(J)=TAUYG(J)/PLAND
800 C**** 6596.
801 C**** ACCUMULATE DIAGNOSTICS 6597.
802 C**** 6598.
803 6000 IM1=I 6662.
804 C**** QUANTITIES ACCUMULATED FOR SURFACE TYPE TABLES IN DIAG1 6663.
805 BLJ(J,9)=BTRHDT
806 BLJ(J,13)=BSHDT
807 BLJ(J,14)=BEVHDT
808 BLJ(J,32)=BTAUL
809 BLJ(J,33)=BTAUF
810 BLJ(J,37)=BWS
811 BLJ(J,28)=BWMG
812 BLJ(J,38)=BTAUL
813 IF(MODDSF.NE.0) GO TO 7000 6673.
814 BLJ(J,23)=BTS
815 BLJ(J,26)=BT2
816 7000 CONTINUE 6677.
817 9000 CONTINUE
818 c write(935),TAU,ELHTG,SHTG,TAUXG,TAUYG
819 C**** 6678.
820 RETURN 6795.
821 990 FORMAT ('0PPBL',3I4,14F8.2) 6818.
822 991 FORMAT ('0SURFACE ',4I4,5F10.4,3F11.7) 6819.
823 992 FORMAT ('0',I2,10F10.4/23X,4F10.4,10X,2F10.4/ 6820.
824 * 33X,3F10.4,10X,2F10.4) 6821.
825 993 FORMAT ('0',I2,10F10.4/23X,7F10.4/33X,7F10.4) 6822.
826 994 FORMAT ('0',I2,11F10.4) 6823.
827 END 6824.

  ViewVC Help
Powered by ViewVC 1.1.22