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

Annotation of /MITgcm_contrib/jscott/igsm/src/surf_land.F

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


Revision 1.1 - (hide annotations) (download)
Fri Aug 11 19:35:32 2006 UTC (18 years, 11 months ago) by jscott
Branch: MAIN
atm2d package

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

  ViewVC Help
Powered by ViewVC 1.1.22