/[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.2 - (hide 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 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 jscott 1.2 #include "CLM.h"
48 jscott 1.1 #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 jscott 1.2 ! common/TSUR/TSURFC(JM0,0:13),TSURFT(JM0),TSURFD(JM0),DTEMSR(JM0)
66     #include "TSRF.COM"
67 jscott 1.1 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 jscott 1.2 WMGMIN=5.
144 jscott 1.1 print *,'WMGMIN 4 LAND=',WMGMIN
145 jscott 1.2 print *,'over land WMG=max(WMG0,WMGMIN)'
146 jscott 1.1 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 jscott 1.2 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 jscott 1.1 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 jscott 1.2 ! WMG=WMG0+WMGMIN
525     ! 07/17/2006
526 jscott 1.1 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 jscott 1.2 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 jscott 1.1 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