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

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

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


Revision 1.3 - (hide annotations) (download)
Mon Apr 23 21:20:18 2007 UTC (18 years, 3 months ago) by jscott
Branch: MAIN
Changes since 1.2: +8 -6 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_OCEAN
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.3 #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     DIMENSION WMGMIN(JM0)
64     LOGICAL POLE,PRNT,HPRNT
65     common/conprn/HPRNT
66 jscott 1.3 ! common/TSUR/TSURFC(JM0,0:13),TSURFT(JM0),TSURFD(JM0),DTEMSR(JM0)
67     #include "TSRF.COM"
68 jscott 1.1 common/SURRAD/TRSURF(JM0,4),SRSURF(JM0,4)
69     c REAL*8 B,TGV,TKV,TSV0,TSV1,TSV 5818.
70     COMMON/CWMG/WMGEA(JM0),NWMGEA(JM0),CHAVER(JM0),DTAV(JM0),DQAV(JM0)
71     & ,Z0AV(JM0),WSAV(JM0),WS0AV(JM0),TAUAV(JM0)
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 jscott 1.3 real T2MZ(JM0)
77 jscott 1.1 C
78     #if ( defined OCEAN_3D || defined ML_2D)
79 jscott 1.2 #include "AGRID.h"
80 jscott 1.1 #endif
81     c
82     DATA RVAP/461.5/ 5819.
83     DATA SHV/0./,SHW/4185./,SHI/2060./,RHOW/1000./,RHOI/916.6/, 5820.
84     * ALAMI/2.1762/,STBO/.5672573E-7/,TF/273.16/,TFO/-1.56/ 5821.
85     DATA Z1I/.1/,Z2LI/2.9/,Z1E/.1/,Z2E/4./,RHOS/91.66/,ALAMS/.35/ 5822.
86     DIMENSION AROUGH(20),BROUGH(20),CROUGH(20),DROUGH(20),EROUGH(20) 5823.
87     DATA AROUGH/16.59,13.99,10.4,7.35,5.241,3.926,3.126,2.632,2.319, 5824.
88     *2.116,1.982,1.893,1.832,1.788,1.757,1.733,1.714,1.699,1.687,1.677/5825.
89     DATA BROUGH/3.245,1.733,0.8481,0.3899,0.1832,0.9026E-1,0.4622E-1, 5826.
90     * .241E-1,.1254E-1,.6414E-2,.3199E-2,.1549E-2,.7275E-3,.3319E-3, 5827.
91     * .1474E-3,.6392E-4,.2713E-4,.1130E-4,.4630E-5,.1868E-5/ 5828.
92     DATA CROUGH/5.111,3.088,1.682,.9239,.5626,.3994,.3282,.3017,.299 5829.
93     *,.3114,.3324,.3587,.3881,.4186,.4492,.4792,.5082,.5361,.5627, 5830.
94     * .5882/ 5831.
95     DATA DROUGH/1.24,1.02,0.806,0.682,0.661,0.771,0.797,0.895,0.994, 5832.
96     * 1.09,1.18,1.27,1.35,1.43,1.50,1.58,1.65,1.71,1.78,1.84/ 5833.
97     DATA EROUGH/0.128,0.130,0.141,0.174,0.238,0.330,0.438,0.550,0.660,5834.
98     * 0.766,0.866,0.962,1.05,1.14,1.22,1.30,1.37,1.45,1.52,1.58/ 5835.
99     QSAT(TM,PR,QLH)=3.797915*EXP(QLH*(7.93252E-6-2.166847E-3/TM))/PR 5836.
100     DLQSDT(TM,QLH)=QLH*2.166847E-3/(TM*TM)
101     c TLOG(Z0)=ALOG(.36*SQRTT/(FMAG*Z0))+2.302585*ROUGH-.08 5837.
102     DATA IFIRST/1/ 5838.
103     ROSNOW(X)=0.54*X/LOG(1.+0.54*X/275.)
104     ALSNOW(X)=2.8E-6*X**2
105     C**** 5839.
106     C**** FDATA 2 LAND COVERAGE (1) 5840.
107     C**** 3 RATIO OF LAND ICE COVERAGE TO LAND COVERAGE (1) 5841.
108     C**** 5842.
109     C**** ODATA 1 OCEAN TEMPERATURE (C) 5843.
110     C**** 2 RATIO OF OCEAN ICE COVERAGE TO WATER COVERAGE (1) 5844.
111     C**** 3 OCEAN ICE AMOUNT OF SECOND LAYER (KG/M**2) 5845.
112     C**** 5846.
113     C**** GDATA 1 OCEAN ICE SNOW AMOUNT (KG/M**2) 5847.
114     C**** 2 EARTH SNOW AMOUNT (KG/M**2) 5848.
115     C**** 3 OCEAN ICE TEMPERATURE OF FIRST LAYER (C) 5849.
116     C**** 4 EARTH TEMPERATURE OF FIRST LAYER (C) 5850.
117     C**** 5 EARTH WATER OF FIRST LAYER (KG/M**2) 5851.
118     C**** 6 EARTH ICE OF FIRST LAYER (KG/M**2) 5852.
119     C**** 7 OCEAN ICE TEMPERATURE OF SECOND LAYER (C) 5853.
120     C**** 8 EARTH TEMPERATURE OF SECOND LAYER (C) 5854.
121     C**** 9 EARTH WATER OF SECOND LAYER (KG/M**2) 5855.
122     C**** 10 EARTH ICE OF SECOND LAYER (KG/M**2) 5856.
123     C**** 12 LAND ICE SNOW AMOUNT (KG/M**2) 5857.
124     C**** 13 LAND ICE TEMPERATURE OF FIRST LAYER (C) 5858.
125     C**** 14 LAND ICE TEMPERATURE OF SECOND LAYER (C) 5859.
126     C**** 5860.
127     C**** BLDATA 1 COMPOSITE SURFACE WIND MAGNITUDE (M/S) 5861.
128     C**** 2 COMPOSITE SURFACE AIR TEMPERATURE (K) 5862.
129     C**** 3 COMPOSITE SURFACE AIR SPECIFIC HUMIDITY (1) 5863.
130     C**** 4 LAYER TO WHICH DRY CONVECTION MIXES (1) 5864.
131     C**** 5 FREE 5865.
132     C**** 6 COMPOSITE SURFACE U WIND 5866.
133     C**** 7 COMPOSITE SURFACE V WIND 5867.
134     C**** 8 COMPOSITE SURFACE MOMENTUM TRANSFER (TAU) 5868.
135     C**** 5869.
136     C**** VDATA 9 WATER FIELD CAPACITY OF FIRST LAYER (KG/M**2) 5870.
137     C**** 10 WATER FIELD CAPACITY OF SECOND LAYER (KG/M**2) 5871.
138     C**** 5872.
139     C**** ROUGHL LOG(ZGS/ROUGHNESS LENGTH) (LOGARITHM TO BASE 10) 5873.
140     C**** ROUGHL will be ROUGHNESS LENGTH
141     C**** 5874.
142     c print *,'surface TAU=',TAU
143     NSTEPS=NSURF*NSTEP/NDYN 5875.
144     IF(IFIRST.NE.1) GO TO 50 5876.
145     print *,' SURFACE after CLM'
146     IFIRST=0 5877.
147     print *,'JM=',jm
148     ! WMGMIN=8.
149     ! WMGMIN=16. ! 2365.05
150     ! WMGMIN=25. ! 2363.05 2366.05
151     WMGM0=8.0
152     WMGM45=25.
153     print *,' WMGM0=', WMGM0,' WMGM45=',WMGM45
154     WMGMAV=0.5*( WMGM0+WMGM45)
155     DWGM=0.5*( WMGM0-WMGM45)
156     do j = 1,jm0
157     rhrad = 3.14159*(-90.+4.*(j-1))/180.
158     WMGMIN(j) = WMGMAV+DWGM*cos(4.*rhrad)
159     enddo
160     print *,'WMGMIN 4 OCEAN is a function of latitude'
161     print 258,(WMGMIN(J),J=1,JM)
162     print *,' WMGE'
163     print 258,(WMGE(1,J),J=1,JM)
164     258 format(12f5.1)
165     COEFSN=100./ROSNOW(10.)
166     COEFSN=1.
167     print *,' COEFSN=',COEFSN
168     do 2567 J=1,JM
169     NWMGEA(J)=0
170     WMGEA(J)=0.
171     CHAVER(J)=0.
172     DTAV(J)=0.
173     DQAV(J)=0.
174     Z0AV(J)=0.
175     WSAV(J)=0.
176     WS0AV(J)=0.
177     TAUAV(J)=0.
178     2567 CONTINUE
179     c10 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     LBLMM1=LBLM-1 5880.
187     IQ1=IM/4+1 5881.
188     IQ2=IM/2+1 5882.
189     IQ3=3*IM/4+1 5883.
190     DTSURF=NDYN*DT/NSURF 5884.
191     print *,' DTSURF=',DTSURF
192     DTSRCE=DT*NDYN 5885.
193     SHA=RGAS/KAPA 5886.
194     RVX=0. 5887.
195     ACE1I=Z1I*RHOI 5888.
196     HC1I=ACE1I*SHI 5889.
197     HC2LI=Z2LI*RHOI*SHI 5890.
198     HC1DE=Z1E*1129950. 5891.
199     HC2DE=Z2E*1129950.+3.5*.125*RHOW*3100. 5892.
200     Z1IBYL=Z1I/ALAMI 5893.
201     Z2LI3L=Z2LI/(3.*ALAMI) 5894.
202     BYRSL=1./(RHOS*ALAMS) 5895.
203     ZS1CO=.5*DSIG(1)*RGAS/GRAV 5896.
204     P1000K=EXPBYK(1000.) 5897.
205     COEFS=GRAV/(100.*DSIG(1)) 5898.
206     COEF1=(1.-SIG(2))/DSIGO(1) 5899.
207     COEF2=(SIG(1)-1.)/DSIGO(1) 5900.
208     DO 20 I=1,IM 5901.
209     SINI(I)=SIN((I-1)*TWOPI/FIM) 5902.
210     20 COSI(I)=COS((I-1)*TWOPI/FIM) 5903.
211     50 S0=S0X*1367./RSDIST 5904.
212     c print *,'IM=',IM
213     c do i=1,200
214     c print *,'HERE1'
215     c enddo
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     TGRND(I,J,2)=GDATA(I,J,3) 5909.
221     c TGRND(I,J,3)=GDATA(I,J,13) 5910.
222     c TGRND(I,J,4)=GDATA(I,J,4) 5911.
223     DO 70 K=1,2 5912.
224     EVAPOR(I,J,K)=0. 5913.
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     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     DO 60 I=1,IM 5924.
239     DU1(I,J)=DUL1(J)
240     60 DV1(I,J)=DVL1(J)
241     C**** 5927.
242     C**** OUTSIDE LOOP OVER J AND I, EXECUTED ONCE FOR EACH GRID POINT 5928.
243     C**** 5929.
244     JPR=-7
245     JPRR=-31
246     DO 7000 J=1,JM 5930.
247     PRNT=j.eq.8
248     PRNT=.FALSE.
249     HEMI=1. 5931.
250     IF(J.LE.JM/2) HEMI=-1. 5932.
251     FCOR=2.*OMEGA*SINP(J) 5933.
252     FMAG=FCOR*HEMI 5934.
253     ROOT2F=SQRT(2.*FMAG) 5935.
254     IF(J.EQ.1) GO TO 80 5936.
255     IF(J.EQ.JM) GO TO 90 5937.
256     WMG0=.5*(WMGE(1,J)+WMGE(1,J+1))+.001 5937.5
257     POLE=.FALSE. 5938.
258     IMAX=IM 5939.
259     GO TO 100 5940.
260     C**** CONDITIONS AT THE SOUTH POLE 5941.
261     80 POLE=.TRUE. 5942.
262     IMAX=1 5943.
263     JVPO=2 5944.
264     RAPO=2.*RAPVN(1) 5945.
265     U1=.25*(U(1,2,1)+V(IQ1,2,1)-U(IQ2,2,1)-V(IQ3,2,1)) 5946.
266     V1=.25*(V(1,2,1)-U(IQ1,2,1)-V(IQ2,2,1)+U(IQ3,2,1)) 5947.
267     WMG0=WMGE(1,2)+.001 5947.5
268     GO TO 100 5948.
269     C**** CONDITIONS AT THE NORTH POLE 5949.
270     90 POLE=.TRUE. 5950.
271     IMAX=1 5951.
272     JVPO=JM 5952.
273     RAPO=2.*RAPVS(JM) 5953.
274     U1=.25*(U(1,JM,1)-V(IQ1,JM,1)-U(IQ2,JM,1)+V(IQ3,JM,1)) 5954.
275     V1=.25*(V(1,JM,1)+U(IQ1,JM,1)-V(IQ2,JM,1)-U(IQ3,JM,1)) 5955.
276     WMG0=WMGE(1,JM)+.001 5955.5
277     C**** ZERO OUT SURFACE DIAGNOSTICS WHICH WILL BE SUMMED OVER LONGITUDE 5956.
278     100 ATRHDT=0. 5957.
279     BTRHDT=0. 5958.
280     CTRHDT=0. 5959.
281     ASHDT=0. 5960.
282     CSHDT=0. 5962.
283     AEVHDT=0. 5963.
284     CEVHDT=0. 5965.
285     ATS=0. 5966.
286     CTS=0. 5968.
287     AT2=0. 5966.
288     CT2=0. 5968.
289     ARIGS=0. 5966.
290     CRIGS=0. 5968.
291     ATAUL=0.
292     ATAUF=0.
293     CTAUL=0.
294     CTAUF=0.
295     AWS=0.
296     CWS=0.
297     AWMG=0.
298     CWMG=0.
299     ACH=0.
300     CCH=0.
301     IM1=IM 5969.
302     DO 6000 I=1,IMAX 5970.
303     C**** 5971.
304     C**** DETERMINE SURFACE CONDITIONS 5972.
305     C**** 5973.
306     PLAND=FDATA(I,J,2) 5974.
307     PWATER=1.-PLAND 5975.
308     PLICE=FDATA(I,J,3)*PLAND 5976.
309     PEARTH=PLAND-PLICE 5977.
310     POICE=ODATA(I,J,2)*PWATER 5978.
311     POCEAN=PWATER-POICE 5979.
312     if(POCEAN.LE.1.E-5)then
313     POCEAN=0.
314     POICE=PWATER
315     endif
316     TTOFR=PEARTH+PLICE+POICE+POCEAN
317     if(abs(TTOFR-1).gt.1.e-3)then
318     print *,' From surface TTOFR=',TTOFR
319     print *,' J=',J,' PLAND=',PLAND,' POCEAN=',POCEAN
320     print *,'POICE=',POICE,' ODATA(I,J,2)=',ODATA(I,J,2)
321     stop
322     end if
323     SP=P(I,J) 5980.
324     PS=SP+PTOP 5981.
325     PSK=EXPBYK(PS) 5982.
326     P1=SIG(1)*SP+PTOP 5983.
327     P1K=EXPBYK(P1) 5984.
328     WSOLD=BLDATA(I,J,1) 5985.
329     USOLD=BLDATA(I,J,6) 5986.
330     VSOLD=BLDATA(I,J,7) 5987.
331     TOLD=BLDATA(I,J,8) 5988.
332     SQRTT=SQRT(TOLD) 5989.
333     GKBYFW=.1296*GRAV/(FCOR*FMAG*WSOLD+1.E-20) 5990.
334     COSWS=GKBYFW*USOLD 5991.
335     SINWS=GKBYFW*VSOLD 5992.
336     IF(POLE) GO TO 1200 5993.
337     U1=.25*(U(IM1,J,1)+U(I,J,1)+U(IM1,J+1,1)+U(I,J+1,1)) 5994.
338     V1=.25*(V(IM1,J,1)+V(I,J,1)+V(IM1,J+1,1)+V(I,J+1,1)) 5995.
339     if(J.eq.JPR.or.J.eq.-12)then
340     print *,' J=',J
341     print *,'TAU=',TAU,' NS=',NS,' ITYPE=',ITYPE
342     print *,'U(I,J,1)=',U(I,J,1),' V(I,J,1)=',V(I,J,1)
343     print *,'U(I,J+1,1)=',U(I,J+1,1),' V(I,J+1,1)=',V(I,J+1,1)
344     print *,'U(IM1,J,1)=',U(IM1,J,1),' V(IM1,J,1)=',V(IM1,J,1)
345     print *,'U(IM1,J+1,1)=',U(IM1,J+1,1),
346     & ' V(IM1,J+1,1)=',V(IM1,J+1,1)
347     endif
348     1200 TH1=T(I,J,1) 5996.
349     Q1=Q(I,J,1) 5997.
350     DTH1=DT1L(J)
351     DQQ1=DQ1L(J)
352     THV1=TH1*(1.+Q1*RVX) 5998.
353     c SRHEAT=SRHR(I,J,1)*COSZ1(I,J)*SRCOR 5999.
354     c SRHDT=SRHEAT*DTSURF 6000.
355     RMBYA=100.*SP*DSIG(1)/GRAV 6001.
356     C**** ZERO OUT QUANTITIES TO BE SUMMED OVER SURFACE TYPES 6002.
357     USS=USSL(J)
358     VSS=VSSL(J)
359     WSS=WSSL(J)
360     TSS=TSSL(J)
361     T2MS=T2ML(J)
362     QSS=QSSL(J)
363     TAUS=TAUSL(J)
364     tauu(j)=tauu(j)+TAUYG(J)*PLAND
365     tauv(j)=tauv(j)+TAUXG(J)*PLAND
366     SINAPS=0. 6009.
367     COSAPS=0. 6010.
368     JR=J
369     DXYPJ=DXYP(J) 6012.
370     TG1S=0. 6013.
371     QGS=0. 6014.
372     BETAS=0. 6015.
373     TRHDTS=0. 6016.
374     SHDTS=0. 6017.
375     EVHDTS=0. 6018.
376     UGS=0. 6019.
377     VGS=0. 6020.
378     WGS=0. 6021.
379     USRS=0. 6022.
380     VSRS=0. 6023.
381     RIS1S=0. 6024.
382     RIGSS=0. 6025.
383     CDMS=0. 6026.
384     CDHS=0. 6027.
385     DGSS=0. 6028.
386     EDS1S=0. 6029.
387     PPBLS=0. 6030.
388     EVAPS=0. 6031.
389     C**** 6032.
390     IF(POCEAN.LE.0.) GO TO 2200 6033.
391     C**** 6034.
392     C**** OCEAN 6035.
393     C**** 6036.
394     ITYPE=1 6037.
395     PTYPE=POCEAN 6038.
396     c formula charnoka
397     TOCEAN=BLDATA(I,J,5)
398     ! ROUGH=MAX(0.018*TOCEAN/GRAV,1.5E-4)
399     ROUGH=MAX(0.035*TOCEAN/GRAV,1.5E-4) ! 2375.05
400     ZGS=10. 6041.
401     NGRNDZ=1 6043.
402     TG1=ODATA(I,J,1) 6044.
403     BETA=1. 6045.
404     ELHX=LHE 6046.
405     TRHT0=TRSURF(J,1)
406     SRHEAT=SRSURF(J,1)*COSZ1(I,J)*SRCOR
407     GO TO 3000 6047.
408     C**** 6048.
409     2200 IF(POICE.LE.0.) GO TO 5000 6049.
410     C**** 6050.
411     C**** OCEAN ICE 6051.
412     C**** 6052.
413     ITYPE=2 6053.
414     PTYPE=POICE 6054.
415     NGRNDZ=NGRND 6055.
416     SNOW=GDATA(I,J,1) 6056.
417     TG1=TGRND(I,J,2) 6057.
418     TG2=GDATA(I,J,7) 6058.
419     ACE2=ODATA(I,J,3) 6059.
420     Z2=ACE2/RHOI 6060.
421     Z2BY4L=Z2/(4.*ALAMI) 6061.
422     if (SNOW.gt.10.)then
423     RHOS0=ROSNOW(SNOW)
424     else
425     RHOS0=275.
426     endif
427     RHOS=COEFSN*RHOS0
428     ALAMS=ALSNOW(RHOS0)
429     BYRSL=1./(RHOS*ALAMS)
430     c Z1BY6L=(Z1IBYL+SNOW*BYRSL)*.1666667 6062.
431     c CDTERM=1.5*TG2-.5*TFO 6063.
432     CDTERM=TG2
433     c CDENOM=1./(2.*Z1BY6L+Z2BY4L) 6064.
434     Z1BY2L=(Z1IBYL+SNOW*BYRSL)*0.5
435     CDENOM=1./(Z1BY2L+2.*Z2BY4L)
436     ROUGH=10**(log10(10.)-4.37)
437     ZGS=10. 6067.
438     HC1=HC1I+SNOW*SHI 6069.
439     BETA=1. 6070.
440     ELHX=LHS 6071.
441     TRHT0=TRSURF(J,3)
442     SRHEAT=SRSURF(J,3)*COSZ1(I,J)*SRCOR
443     GO TO 3000 6072.
444     C**** 6073.
445     C
446     C**** BOUNDARY LAYER INTERACTION 6135.
447     C**** 6136.
448     3000 continue
449    
450     SRHDT=SRHEAT*DTSURF
451     TKV=THV1*PSK 6137.
452     ZS1=ZS1CO*TKV*SP/PS 6138.
453     P1=SIG(1)*SP+PTOP 6139.
454     DTGRND=DTSURF/NGRNDZ 6143.
455     SHDT=0. 6144.
456     EVHDT=0. 6145.
457     TRHDT=0. 6146.
458     F1DT=0. 6147.
459     C**** LOOP OVER GROUND TIME STEPS 6148.
460     c DO 3600 NG=1,NGRNDZ 6149.
461     TG=TG1+TF 6150.
462     QG=QSAT(TG,PS,ELHX) 6151.
463     TGV=TG*(1.+QG*RVX) 6152.
464     UG=U1
465     VG=V1
466     UG=0.75*U1 ! 2372
467     VG=0.75*V1 ! 2372
468     W1=SQRT(U1*U1+V1*V1)
469     W1=SQRT(UG*UG+VG*VG) ! 2372
470     WS0=W1
471     WMG=WMG0+WMGMIN(J)
472     WS=SQRT((0.75*W1)**2+WMG)
473     WS=SQRT(W1**2+WMG) ! 2372
474     if(J.eq.JPR)then
475     print *,' '
476     print *,'TAU=',TAU,' NS=',NS,' ITYPE=',ITYPE
477     print *,'TG=',TG,' QG=',QG
478     print *,'RVX=',RVX,' TG1=',TG1
479     endif
480    
481     #if ( defined CPL_OCEANCO2 || defined OCEAN_3D)
482     if(ITYPE.eq.1)then
483     NWMGEA(J)=NWMGEA(J)+1
484     WSAV(J)=WSAV(J)+WS
485     end if
486     #endif
487    
488     ! WG=WS 2369
489     WG=WS
490     THS=TH1
491     QS=Q1
492     TSV=THS*PSK
493     Z0=ROUGH
494     ROUGH=log10(ZGS/ROUGH)
495     CDN=.0231/(ROUGH*ROUGH)
496     LR=ROUGH*2.-.5
497     IF(LR.GT.20) LR=20
498     IF(LR.LT.1) LR=1
499     RIGS=ZGS*GRAV*(TSV-TGV)/(TGV*WS*WS)
500     SINAP=0.
501     COSAP=1.
502     IF(RIGS.LE.0) THEN
503     C surface layer has unstable stratification
504     CIA=TWOPI*0.0625/(1.+WS*CIAX)
505     DM=SQRT((1.-AROUGH(LR)*RIGS)*(1.-BROUGH(LR)*RIGS)/
506     * (1.-CROUGH(LR)*RIGS))
507     DH=1.35*SQRT((1.-DROUGH(LR)*RIGS)/(1.-EROUGH(LR)*RIGS))
508     ELSE
509     C surface layer has stable stratification
510     CIA=TWOPI*(0.09375-0.03125/(1.+4*RIGS**2))/(1.+WS*CIAX)
511     DM=1./(1.+(11.238+89.9*RIGS)*RIGS)
512     DH=1.35/(1.+1.93*RIGS)
513     END IF
514     CDH=CDN*DM*DH
515     USR=COS(CIA)
516     VSR=SIN(CIA)*HEMI
517     US=(USR*UG-VSR*VG)
518     VS=(VSR*UG+USR*VG)
519     if(J.eq.JPRR)then
520     print *,'TAU=',TAU,' NS=',NS,' ITYPE=',ITYPE
521     print *,'WS=',WS,' ZGS=',ZGS
522     print *,'DM=',DM,' DH=',DH
523     print *,'RIGS=',RIGS,' TGV=',TGV
524     endif
525     RCDHWS=CDH*WS*100.*PS/(RGAS*TSV)
526     if(J.eq.JPR)then
527     c print *,' '
528     print *,'TAU=',TAU,' NS=',NS,' ITYPE=',ITYPE
529     print *,'CDH=',CDH,' RGAS=',RGAS
530     print *,'PS=',PS,' TSV=',TSV
531     print *,'WS=',WS,' RCDHWS=',RCDHWS
532     endif
533     TS=TSV/(1.+QS*RVX) 6467.
534     QSATS=QSAT(TS,PS,ELHX) 6468.
535     c dLQS/dTs
536     DLQSDTS=DLQSDT(TS,ELHX)
537     c dLQS/dTs
538     IF(QS.LE.QSATS) GO TO 3500 6469.
539     DQSSDT=QSATS*ELHX/(RVAP*TS*TS) 6470.
540     X=(QS-QSATS)/(DQSSDT+(SHA/ELHX)) 6471.
541     TS=TS+X 6472.
542     QS=QS+X*(SHA/ELHX) 6473.
543     C**** CALCULATE RHOS*CDM*WS AND RHOS*CDH*WS 6474.
544     3500 CDM=CDN*DM 6475.
545     RCDMWS=CDM*WS*100.*PS/(RGAS*TS) 6476.
546     C**** CALCULATE FLUXES OF SENSIBLE HEAT, LATENT HEAT, THERMAL 6478.
547     C**** RADIATION, AND CONDUCTION HEAT (WATTS/M**2) 6479.
548     SHEAT=SHA*RCDHWS*(TS-TG) 6480.
549     BETAUP=BETA 6481.
550     IF(QS.GT.QG) BETAUP=1. 6482.
551     EVHEAT=(LHE+TG1*SHV)*BETAUP*RCDHWS*(QS-QG) 6483.
552     TRHEAT=TRHT0-STBO*(TG*TG)*(TG*TG)
553     if(J.eq.JPRR)then
554     c print *,' '
555     ! print *,'TAU=',TAU,' NS=',NS,' ITYPE=',ITYPE
556     ! print *,'TRHT0=',TRHT0,' STBO=',STBO
557     ! print *,'TG=',TG,' TS=',TS
558     ! print *,'TRHEAT=',TRHEAT
559     print *,'EVHEAT=',EVHEAT
560     ! print *,'SHA=',SHA,' RCDHWS=',RCDHWS
561     print *,'SHEAT=',SHEAT
562     endif
563     TG1OLD=TG1
564     SHEATOLD=SHEAT
565     #if ( defined OCEAN_3D )
566     IF(ITYPE.EQ.1 .or. ITYPE.EQ.2) GO TO 3620
567     #else
568     IF(ITYPE.EQ.1) GO TO 3620 6485.
569     #endif
570     C**** CALCULATE FLUXES USING IMPLICIT TIME STEP FOR NON-OCEAN POINTS 6486.
571     F0=SRHEAT+TRHEAT+SHEAT+EVHEAT 6487.
572     c F1=(TG1-CDTERM-F0*Z1BY6L)*CDENOM 6488.
573     F1=(TG1-CDTERM)*CDENOM
574     DSHDTG=-RCDHWS*SHA
575     DQGDTG=QG*ELHX/(RVAP*TG*TG) 6490.
576     DEVDTG=-RCDHWS*LHE*BETAUP*DQGDTG
577     DTRDTG=-4.*STBO*TG*TG*TG 6492.
578     DF0DTG=DSHDTG+DEVDTG+DTRDTG 6493.
579     c DFDTG=DF0DTG-(1.-DF0DTG*Z1BY6L)*CDENOM 6493.5
580     DFDTG=DF0DTG-CDENOM
581     c DF1DTG=(1.-DF0DTG*Z1BY6L)*CDENOM
582     DF1DTG=CDENOM
583     DTG=(F0-F1)*DTGRND/(HC1-DTGRND*DFDTG) 6494.
584     SHDT=SHDT+DTGRND*(SHEAT+DTG*DSHDTG) 6495.
585     EVHDT=EVHDT+DTGRND*(EVHEAT+DTG*DEVDTG) 6496.
586     TRHDT=TRHDT+DTGRND*(TRHEAT+DTG*DTRDTG) 6497.
587     TG1=TG1+DTG
588     c F1DT=F1DT+DTGRND*(TG1-CDTERM-(F0+DTG*DF0DTG)*Z1BY6L)*CDENOM 6498.
589     F1DT=F1DT+DTGRND*(TG1-CDTERM)*CDENOM
590     DU1(I,J)=DU1(I,J)+PTYPE*DTGRND*RCDMWS*US*COEFS/SP 6499.
591     DV1(I,J)=DV1(I,J)+PTYPE*DTGRND*RCDMWS*VS*COEFS/SP 6500.
592     c TG1=TG1+DTG 6501.
593     3600 CONTINUE 6502.
594     GO TO 3700 6503.
595     C**** CALCULATE FLUXES USING EXPLICIT TIME STEP FOR OCEAN POINTS 6504.
596     3620 SHDT=DTSURF*SHEAT 6505.
597     EVHDT=DTSURF*EVHEAT 6506.
598     TRHDT=DTSURF*TRHEAT 6507.
599     DU1(I,J)=DU1(I,J)+PTYPE*DTSURF*RCDMWS*US*COEFS/SP 6508.
600     DV1(I,J)=DV1(I,J)+PTYPE*DTSURF*RCDMWS*VS*COEFS/SP 6509.
601     3700 CONTINUE
602     EPS=1.D-8
603     c print *,'FROM SURFACE NS=',NS
604     c print *,'J=',J,' ITYPE=',ITYPE
605     c print *,RCDMWS,WS
606     WWS=max(W1,1.D-4)
607     ZTEM=ZGS
608     ZTEM=2.0
609     CALL TZM(T2M,TH2M,ZTEM,Z0,ZGS,SP,TG,TS,RIGS,WS,
610     & -SHEATOLD,RCDMWS*WS,LR,EPS)
611     c print *,'FROM SURFACE'
612     c print *,'J=',J,' ITYPE=',ITYPE
613     c print *,'POCEAN=',POCEAN,' PTYPE=',PTYPE
614     c print *,'TS=',TS,' TG=',TG
615     c print *,' T2M=',T2M,' TH2M=',TH2M
616     F0DT=CORSR*SRHDT+TRHDT+SHDT+EVHDT 6510.
617     if(J.eq.JPR)then
618     print *,'TAU=',TAU,' NS=',NS,' ITYPE=',ITYPE
619     print *,'DTSURF=',DTSURF,' CORSR=',CORSR
620     print *,'SRHDT=',SRHDT,' TRHDT=',TRHDT
621     print *,'SHDT=',SHDT,' EVHDT=',EVHDT
622     print *,'F0DT=',F0DT
623     print *,'US=',US,' VS=',VS
624     print *,'COEFS=',COEFS,' SP=',SP
625     endif
626     c print *,'From surface ',TAU,CORSR,SRHDT,TRHDT,SHDT,EVHDT
627     C**** CALCULATE EVAPORATION 6511.
628     CCC DQ1=EVHDT/((LHE+TG1*SHV)*RMBYA) 6512.
629     DQ1=EVHDT/(ELHX*RMBYA)
630     IF(DQ1*PTYPE.LE.Q1) GO TO 3720 6513.
631     DQ1=Q1/PTYPE 6514.
632     CCC EVHDT=DQ1*(LHE+TG1*SHV)*RMBYA 6515.
633     EVHDT=DQ1*ELHX*RMBYA
634     3720 EVAP=-DQ1*RMBYA 6516.
635     C**** ACCUMULATE SURFACE FLUXES AND PROGNOSTIC AND DIAGNOSTIC QUANTITIES6517.
636     E0(I,J,ITYPE)=E0(I,J,ITYPE)+F0DT 6518.
637     E1(I,J,ITYPE)=E1(I,J,ITYPE)+F1DT 6519.
638     EVAPOR(I,J,ITYPE)=EVAPOR(I,J,ITYPE)+EVAP 6520.
639     if(PRNT)then
640     c write(78,*) ,' '
641     c write(78,*) ,'TAU=',TAU
642     write(78,*) ,'J=',j,' ITYPE=',ITYPE,' PTYPE=',PTYPE
643     write(78,*) ,'TS=',TS,' TG=',TG,' QS=',QS
644     write(78,*) ,'TG1=',TG1,' TG1OLD=',TG1OLD
645     write(78,*) ,'TG2=',TG2
646     write(78,*) ,'SHEAT=',SHEAT,' EVHEAT=',EVHEAT
647     write(78,*) ,'TRHEAT=',TRHEAT,' SRHEAT=',SRHEAT
648     write(78,*) ,'EVAP mm/day=',24.*3600.*EVAP/DTSURF
649     write(78,*) ,'EVAP=',EVAP,
650     & ' F0DT=',F0DT/DTSURF,' F1DT=',F1DT/DTSURF
651     endif
652     #if ( defined OCEAN_3D || defined ML_2D )
653     C For ocean model
654     #if ( defined ML_2D )
655     if(ITYPE.eq.1)then
656     #endif
657     C DNetHeat by DTG
658     DSHDTG=-RCDHWS*SHA
659     DQGDTG=QG*ELHX/(RVAP*TG*TG)
660     DEVDTG=-RCDHWS*LHE*BETAUP*DQGDTG
661     DTRDTG=-4.*STBO*TG*TG*TG
662     DF0DTG=DSHDTG+DEVDTG+DTRDTG
663     if(EVHEAT.lt.0.0)then
664     DEVDTGEQ=EVHEAT*DLQSDTS
665     else
666     DEVDTGEQ=0.0
667     endif
668     C DNetHeat by DTG
669     #if ( defined OCEAN_3D )
670     if(ITYPE.eq.1)then
671     #endif
672     dhfodtg(j)=dhfodtg(j)+DF0DTG
673     devodtg(j)=devodtg(j)-DEVDTG/LHE
674     dhfodtgeq(j)=dhfodtgeq(j)+DEVDTGEQ
675     devodtgeq(j)=devodtgeq(j)-DEVDTGEQ/LHE
676     evao(j)=evao(j)+EVAP
677     hfluxo(j)=hfluxo(j)+F0DT
678     naveo(j)=naveo(j)+1
679     endif
680     if(ITYPE.eq.2)then
681     evai(j)=evai(j)+EVAP
682     hfluxi(j)=hfluxi(j)+F0DT
683     dhfidtg(j)=dhfidtg(j)+DF0DTG
684     devidtg(j)=devidtg(j)-DEVDTG/LHE
685     dhfidtgeq(j)=dhfidtgeq(j)+DEVDTGEQ
686     devidtgeq(j)=devidtgeq(j)-DEVDTGEQ/LHE
687     c tairi(j)=tairi(j)+TS
688     navei(j)=navei(j)+1
689     endif
690     tauu(j)=tauu(j)+RCDMWS*US*PTYPE
691     tauv(j)=tauv(j)+RCDMWS*VS*PTYPE
692     C For ocean model
693     #endif
694     TGRND(I,J,ITYPE)=TG1 6521.
695     TSSFC(I,J,ITYPE)=TS 6521.5
696     c TH1=TH1-SHDT*PTYPE/(SHA*RMBYA*P1K) 6522.
697     c Q1=Q1-DQ1*PTYPE 6523.
698    
699     DTH1=DTH1-SHDT*PTYPE/(SHA*RMBYA*P1K)
700     DQQ1=DQQ1-DQ1*PTYPE
701    
702     USS=USS+US*PTYPE 6524.
703     VSS=VSS+VS*PTYPE 6525.
704     WSS=WSS+WS*PTYPE 6526.
705     TSS=TSS+TS*PTYPE 6527.
706     T2MS=T2MS+T2M*PTYPE
707     QSS=QSS+QS*PTYPE 6528.
708     ! TAUS=TAUS+CDM*WS*W1*PTYPE 6529.
709     TAUS=TAUS+CDM*WS*WS*PTYPE ! 2373
710     SINAPS=SINAPS+SINAP*PTYPE 6530.
711     COSAPS=COSAPS+COSAP*PTYPE 6531.
712     TG1S=TG1S+TG1*PTYPE 6532.
713     QGS=QGS+QG*PTYPE 6533.
714     BETAS=BETAS+BETA*PTYPE 6534.
715     TRHDTS=TRHDTS+TRHDT*PTYPE 6535.
716     SHDTS=SHDTS+SHDT*PTYPE 6536.
717     EVHDTS=EVHDTS+EVHDT*PTYPE 6537.
718     UGS=UGS+UG*PTYPE 6538.
719     VGS=VGS+VG*PTYPE 6539.
720     WGS=WGS+WG*PTYPE 6540.
721     USRS=USRS+USR*PTYPE 6541.
722     VSRS=VSRS+VSR*PTYPE 6542.
723     c RIS1S=RIS1S+RIS1*PTYPE 6543.
724     RIGSS=RIGSS+RIGS*PTYPE 6544.
725     CDMS=CDMS+CDM*PTYPE 6545.
726     CDHS=CDHS+CDH*PTYPE 6546.
727     c DGSS=DGSS+DGS*PTYPE 6547.
728     c EDS1S=EDS1S+EDS1*PTYPE 6548.
729     c PPBLS=PPBLS+PPBL*PTYPE 6549.
730     EVAPS=EVAPS+EVAP*PTYPE 6550.
731     GO TO (4000,4100,5000,5000),ITYPE 6551.
732     C**** 6552.
733     C**** OCEAN 6553.
734     C**** 6554.
735     4000 ASHDT=ASHDT+SHDT*POCEAN 6555.
736     AEVHDT=AEVHDT+EVHDT*POCEAN 6556.
737     ATRHDT=ATRHDT+TRHDT*POCEAN 6557.
738     ATS=ATS+(TS-TF)*POCEAN 6558.
739     c AT2=AT2+(TH2M-TF)*POCEAN
740     AT2=AT2+(T2M-TF)*POCEAN
741     ARIGS=ARIGS+RIGS*POCEAN
742     BLDATA(I,J,5)=CDM*WS*W1
743     BLDATA(I,J,5)=CDM*WS*WS ! 2373
744     ATAUL=ATAUL+RCDMWS*US*POCEAN
745     ATAUF=ATAUF+RCDMWS*VS*POCEAN
746     AWS=AWS+WS*POCEAN
747     AWMG=AWMG+SQRT(WMG)*POCEAN
748     ACH=ACH+RCDHWS*POCEAN
749     GO TO 2200 6559.
750     C**** 6560.
751     C**** OCEAN ICE 6561.
752     C**** 6562.
753     4100 CSHDT=CSHDT+SHDT*POICE 6563.
754     CEVHDT=CEVHDT+EVHDT*POICE 6564.
755     CTRHDT=CTRHDT+TRHDT*POICE 6565.
756     CTS=CTS+(TS-TF)*POICE 6566.
757     c CT2=CT2+(TH2M-TF)*POICE 6566.
758     CT2=CT2+(T2M-TF)*POICE
759     CRIGS=CRIGS+RIGS*POICE
760     CTAUL=CTAUL+RCDMWS*US*POICE
761     CTAUF=CTAUF+RCDMWS*VS*POICE
762     CWS=CWS+WS*POICE
763     CWMG=CWMG+SQRT(WMG)*POICE
764     CCH=CCH+RCDHWS*POICE
765     GO TO 5000 6567.
766     C**** 6568.
767     C**** NON-OCEAN POINTS WHICH ARE NOT MELTING OR FREEZING WATER USE 6583.
768     C**** IMPLICIT TIME STEPS 6584.
769     C**** 6585.
770     C**** UPDATE SURFACE AND FIRST LAYER QUANTITIES 6586.
771     C**** 6587.
772     5000 CONTINUE
773     c print *,j,TH1,Q1
774     c print *,j,DU1(1,j),DV1(1,j)
775     T(I,J,1)=TH1 6588.
776     & +DTH1
777     Q(I,J,1)=Q1 6589.
778     & +DQQ1
779     BLDATA(I,J,1)=WSS 6590.
780     BLDATA(I,J,2)=TSS 6591.
781     BLDATA(I,J,3)=QSS 6592.
782     BLDATA(I,J,6)=USS 6593.
783     BLDATA(I,J,7)=VSS 6594.
784     BLDATA(I,J,8)=TAUS 6595.
785 jscott 1.3 T2MZ(J)=T2MS
786 jscott 1.1 c print *,j,T(I,J,1),Q(I,J,1)
787     c print *,(TGRND(I,J,k),k=1,4)
788     c print *,(EVAPOR(I,J,k),k=1,4)
789     c print *,(E0(I,J,k),k=1,4)
790     c print *,(E1(I,J,k),k=1,4)
791     c print *,j,DU1(1,j),DV1(1,j)
792     C**** 6596.
793     C**** ACCUMULATE DIAGNOSTICS 6597.
794     C**** 6598.
795     C**** QUANTITIES ACCUMULATED FOR REGIONS IN DIAG1 6599.
796     IF(JR.EQ.JM) GO TO 5700 6600.
797     DJ(JR,9)=DJ(JR,9)+TRHDTS*DXYPJ 6601.
798     DJ(JR,13)=DJ(JR,13)+SHDTS*DXYPJ 6602.
799     DJ(JR,14)=DJ(JR,14)+EVHDTS*DXYPJ 6603.
800     DJ(JR,19)=DJ(JR,19)+EVAPS*DXYPJ 6604.
801     IF(MODDSF.NE.0) GO TO 5700 6605.
802     DJ(JR,23)=DJ(JR,23)+(TSS-TF)*DXYPJ 6606.
803     5700 CONTINUE
804     6000 IM1=I 6662.
805     C**** QUANTITIES ACCUMULATED FOR SURFACE TYPE TABLES IN DIAG1 6663.
806     AJ(J,9)=AJ(J,9)+ATRHDT 6664.
807     BJ(J,9)=BJ(J,9)+BLJ(J,9)
808     CJ(J,9)=CJ(J,9)+CTRHDT 6666.
809     AJ(J,13)=AJ(J,13)+ASHDT 6667.
810     BJ(J,13)=BJ(J,13)+BLJ(J,13)
811     CJ(J,13)=CJ(J,13)+CSHDT 6669.
812     AJ(J,14)=AJ(J,14)+AEVHDT 6670.
813     BJ(J,14)=BJ(J,14)+BLJ(J,14)
814     if(J.eq.-23)then
815     print *,'TAU=',TAU,' BJ(J,14)=',BJ(J,14)
816     endif
817     CJ(J,14)=CJ(J,14)+CEVHDT 6672.
818     AJ(J,32)=AJ(J,32)+ATAUL
819     BJ(J,32)=BJ(J,32)+BLJ(J,32)
820     CJ(J,32)=CJ(J,32)+CTAUL
821     AJ(J,33)=AJ(J,33)+ATAUF
822     BJ(J,33)=BJ(J,33)+BLJ(J,33)
823     CJ(J,33)=CJ(J,33)+CTAUF
824     AJ(J,37)=AJ(J,37)+AWS
825     BJ(J,37)=BJ(J,37)+BLJ(J,37)
826     CJ(J,37)=CJ(J,37)+CWS
827     AJ(J,28)=AJ(J,28)+AWMG
828     BJ(J,28)=BJ(J,28)+BLJ(J,28)
829     CJ(J,28)=CJ(J,28)+CWMG
830     AJ(J,38)=AJ(J,38)+ATAUL
831     BJ(J,38)=BJ(J,38)+BLJ(J,38)
832     CJ(J,38)=CJ(J,38)+CTAUL
833     IF(MODDSF.NE.0) GO TO 7000 6673.
834     AJ(J,23)=AJ(J,23)+ATS 6674.
835     BJ(J,23)=BJ(J,23)+BLJ(J,23)
836     CJ(J,23)=CJ(J,23)+CTS 6676.
837     AJ(J,26)=AJ(J,26)+AT2 6674.
838     BJ(J,26)=BJ(J,26)+BLJ(J,26)
839     CJ(J,26)=CJ(J,26)+CT2 6676.
840     AJ(J,27)=AJ(J,27)+ARIGS 6674.
841     CJ(J,27)=CJ(J,27)+CRIGS 6676.
842    
843     c print *,j,'CTS=',CTS,' CT2=',CT2
844     c print *,'BLDATA'
845     c print *,(BLDATA(1,j,k),k=1,3)
846     c print *,(BLDATA(1,j,k),k=6,8)
847     7000 CONTINUE 6677.
848     C**** 6678.
849     C**** ADD IN SURFACE FRICTION TO FIRST LAYER WIND 6679.
850     C**** 6680.
851     DO 7600 I=1,IM 6681.
852     U(I,2,1)=U(I,2,1)-2.*(DU1(1,1)*COSI(I)-DV1(1,1)*SINI(I))*RAPVN(1) 6682.
853     V(I,2,1)=V(I,2,1)-2.*(DV1(1,1)*COSI(I)+DU1(1,1)*SINI(I))*RAPVN(1) 6683.
854     U(I,JM,1)=U(I,JM,1) 6684.
855     * -2.*(DU1(1,JM)*COSI(I)+DV1(1,JM)*SINI(I))*RAPVS(JM) 6685.
856     7600 V(I,JM,1)=V(I,JM,1) 6686.
857     * -2.*(DV1(1,JM)*COSI(I)-DU1(1,JM)*SINI(I))*RAPVS(JM) 6687.
858     DO 7700 J=2,JMM1 6688.
859     I=IM 6689.
860     DO 7700 IP1=1,IM 6690.
861     if(J.eq.JPR.or.J.eq.-12)then
862     print *,' J=',J,' before'
863     print *,'U(I,J,1)=',U(I,J,1),' V(I,J,1)=',V(I,J,1)
864     print *,'U(I,J+1,1)=',U(I,J+1,1),' V(I,J+1,1)=',V(I,J+1,1)
865     print *,'DU1(I,J)=',DU1(I,J),' DU1(IP1,J)=',DU1(IP1,J)
866     endif
867     U(I,J,1)=U(I,J,1)-(DU1(I,J)+DU1(IP1,J))*RAPVS(J) 6691.
868     V(I,J,1)=V(I,J,1)-(DV1(I,J)+DV1(IP1,J))*RAPVS(J) 6692.
869     U(I,J+1,1)=U(I,J+1,1)-(DU1(I,J)+DU1(IP1,J))*RAPVN(J) 6693.
870     V(I,J+1,1)=V(I,J+1,1)-(DV1(I,J)+DV1(IP1,J))*RAPVN(J) 6694.
871     if(J.eq.JPR.or.J.eq.-12)then
872     print *,' J=',J,' after'
873     print *,'U(I,J,1)=',U(I,J,1),' V(I,J,1)=',V(I,J,1)
874     print *,'U(I,J+1,1)=',U(I,J+1,1),' V(I,J+1,1)=',V(I,J+1,1)
875     print *,'DU1(I,J)=',DU1(I,J),' DU1(IP1,J)=',DU1(IP1,J)
876     endif
877     7700 I=IP1 6695.
878     c print *,'U V'
879     c do j=1,jm
880     c print *,j,U(I,J,1),v(I,J,1)
881     c enddo
882     C**** 6696.
883     C**** DRY CONVECTION ORIGINATING FROM THE FIRST LAYER 6697.
884     C**** 6698.
885     C**** LOAD U,V INTO UT,VT. UT,VT WILL BE FIXED DURING DRY CONVECTION 6699.
886     C**** WHILE U,V WILL BE UPDATED. 6700.
887     DO 8050 L=1,LM 6701.
888     DO 8050 J=2,JM 6702.
889     DO 8050 I=1,IM 6703.
890     UT(I,J,L)=U(I,J,L) 6704.
891     8050 VT(I,J,L)=V(I,J,L) 6705.
892     C**** OUTSIDE LOOPS OVER J AND I 6706.
893     DO 8500 J=1,JM 6707.
894     POLE=.FALSE. 6708.
895     IF(J.EQ.1.OR.J.EQ.JM) POLE=.TRUE. 6709.
896     IMAX=IM 6710.
897     IF(POLE) IMAX=IM 6711.
898     DO 8120 K=1,2 6712.
899     RA(K)=RAPVS(J) 6713.
900     8120 RA(K+2)=RAPVN(J) 6714.
901     IM1=IM 6715.
902     DO 8500 I=1,IMAX 6716.
903     BLDATA(I,J,4)=1. 6717.
904     IF(T(I,J,1)*(1.+Q(I,J,1)*RVX).LE. 6718.
905     * T(I,J,2)*(1.+Q(I,J,2)*RVX)) GO TO 8500 6719.
906     C**** MIX HEAT AND MOISTURE THROUGHOUT THE BOUNDARY LAYER 6720.
907     PKMS=PK(I,J,1)*DSIG(1)+PK(I,J,2)*DSIG(2) 6721.
908     THPKMS=T(I,J,1)*(PK(I,J,1)*DSIG(1))+T(I,J,2)*(PK(I,J,2)*DSIG(2)) 6722.
909     QMS=Q(I,J,1)*DSIG(1)+Q(I,J,2)*DSIG(2) 6723.
910     TVMS=T(I,J,1)*(1.+Q(I,J,1)*RVX)*(PK(I,J,1)*DSIG(1)) 6724.
911     * +T(I,J,2)*(1.+Q(I,J,2)*RVX)*(PK(I,J,2)*DSIG(2)) 6725.
912     THETA=TVMS/PKMS 6726.
913    
914     #if ( defined CPL_CHEM )
915     !
916     ! --- 03/23/95
917     !
918     cfc11ms = cfc11(i,j,1)*dsig(1) + cfc11(i,j,2)*dsig(2)
919    
920     cfc12ms = cfc12(i,j,1)*dsig(1) + cfc12(i,j,2)*dsig(2)
921    
922     xn2oms = xn2o (i,j,1)*dsig(1) + xn2o (i,j,2)*dsig(2)
923    
924     o3ms = o3 (i,j,1)*dsig(1) + o3 (i,j,2)*dsig(2)
925    
926     coms = co (i,j,1)*dsig(1) + co (i,j,2)*dsig(2)
927    
928     zco2ms = zco2 (i,j,1)*dsig(1) + zco2 (i,j,2)*dsig(2)
929    
930     xnoms = xno (i,j,1)*dsig(1) + xno (i,j,2)*dsig(2)
931    
932     xno2ms = xno2 (i,j,1)*dsig(1) + xno2 (i,j,2)*dsig(2)
933    
934     xn2o5ms = xn2o5(i,j,1)*dsig(1) + xn2o5(i,j,2)*dsig(2)
935    
936     hno3ms = hno3 (i,j,1)*dsig(1) + hno3 (i,j,2)*dsig(2)
937    
938     ch4ms = ch4 (i,j,1)*dsig(1) + ch4 (i,j,2)*dsig(2)
939    
940     ch2oms = ch2o (i,j,1)*dsig(1) + ch2o (i,j,2)*dsig(2)
941    
942     so2ms = so2 (i,j,1)*dsig(1) + so2 (i,j,2)*dsig(2)
943    
944     h2so4ms = h2so4(i,j,1)*dsig(1) + h2so4(i,j,2)*dsig(2)
945    
946     ! === if hfc, pfc, and sf6 are included:
947     #ifdef INC_3GASES
948     ! === 032698
949     hfc134ams = hfc134a(i,j,1)*dsig(1)
950     & + hfc134a(i,j,2)*dsig(2)
951    
952     pfcms = pfc(i,j,1)*dsig(1)
953     & + pfc(i,j,2)*dsig(2)
954    
955     sf6ms = sf6(i,j,1)*dsig(1)
956     & + sf6(i,j,2)*dsig(2)
957     ! ===
958     #endif
959    
960     bcms = bcarbon (i,j,1)*dsig(1) + bcarbon (i,j,2)*dsig(2)
961     ocms = ocarbon (i,j,1)*dsig(1) + ocarbon (i,j,2)*dsig(2)
962    
963     c 062295
964     c h2o2ms = h2o2 (i,j,1)*dsig(1) + h2o2 (i,j,2)*dsig(2)
965    
966     !
967     #endif
968    
969     DO 8140 L=3,LM 6727.
970     IF(THETA.LT.T(I,J,L)*(1.+Q(I,J,L)*RVX)) GO TO 8160 6728.
971     PKMS=PKMS+(PK(I,J,L)*DSIG(L)) 6729.
972     THPKMS=THPKMS+T(I,J,L)*(PK(I,J,L)*DSIG(L)) 6730.
973     QMS=QMS+Q(I,J,L)*DSIG(L) 6731.
974     TVMS=TVMS+T(I,J,L)*(1.+Q(I,J,L)*RVX)*(PK(I,J,L)*DSIG(L)) 6732.
975    
976     #if ( defined CPL_CHEM )
977     !
978     ! --- 03/23/95
979     !
980     cfc11ms = cfc11ms + cfc11(i,j,l)*dsig(l)
981    
982     cfc12ms = cfc12ms + cfc12(i,j,l)*dsig(l)
983    
984     xn2oms = xn2oms + xn2o (i,j,l)*dsig(l)
985    
986     o3ms = o3ms + o3 (i,j,l)*dsig(l)
987    
988     coms = coms + co (i,j,l)*dsig(l)
989    
990     zco2ms = zco2ms + zco2 (i,j,l)*dsig(l)
991    
992     xnoms = xnoms + xno (i,j,l)*dsig(l)
993    
994     xno2ms = xno2ms + xno2 (i,j,l)*dsig(l)
995    
996     xn2o5ms = xn2o5ms + xn2o5(i,j,l)*dsig(l)
997    
998     hno3ms = hno3ms + hno3 (i,j,l)*dsig(l)
999    
1000     ch4ms = ch4ms + ch4 (i,j,l)*dsig(l)
1001    
1002     ch2oms = ch2oms + ch2o (i,j,l)*dsig(l)
1003    
1004     so2ms = so2ms + so2 (i,j,l)*dsig(l)
1005    
1006     h2so4ms = h2so4ms + h2so4(i,j,l)*dsig(l)
1007    
1008     ! === if hfc, pfc, and sf6 are included:
1009     #ifdef INC_3GASES
1010     ! === 032698
1011     hfc134ams = hfc134ams
1012     & + hfc134a(i,j,l)*dsig(l)
1013    
1014     pfcms = pfcms
1015     & + pfc(i,j,l)*dsig(l)
1016    
1017     sf6ms = sf6ms
1018     & + sf6(i,j,l)*dsig(l)
1019     ! ===
1020     #endif
1021    
1022     bcms = bcms + bcarbon (i,j,l)*dsig(l)
1023     ocms = ocms + ocarbon (i,j,l)*dsig(l)
1024    
1025     c 062295
1026     c h2o2ms = h2o2ms + h2o2 (i,j,l)*dsig(l)
1027     !
1028     #endif
1029    
1030     8140 THETA=TVMS/PKMS 6733.
1031     L=LM+1 6734.
1032     8160 LMAX=L-1 6735.
1033     RDSIGS=1./(SIGE(1)-SIGE(LMAX+1)) 6736.
1034     THM=THPKMS/PKMS 6737.
1035     QMS=QMS*RDSIGS 6738.
1036    
1037     #if ( defined CPL_CHEM )
1038     !
1039     ! --- 03/23/95
1040     !
1041     cfc11ms = cfc11ms*rdsigs
1042    
1043     cfc12ms = cfc12ms*rdsigs
1044    
1045     xn2oms = xn2oms *rdsigs
1046    
1047     o3ms = o3ms *rdsigs
1048    
1049     coms = coms *rdsigs
1050    
1051     zco2ms = zco2ms *rdsigs
1052    
1053     xnoms = xnoms *rdsigs
1054    
1055     xno2ms = xno2ms *rdsigs
1056    
1057     xn2o5ms = xn2o5ms*rdsigs
1058    
1059     hno3ms = hno3ms *rdsigs
1060    
1061     ch4ms = ch4ms *rdsigs
1062    
1063     ch2oms = ch2oms *rdsigs
1064    
1065     so2ms = so2ms *rdsigs
1066    
1067     h2so4ms = h2so4ms*rdsigs
1068    
1069     ! === if hfc, pfc, and sf6 are included:
1070     #ifdef INC_3GASES
1071     ! === 032698
1072     hfc134ams = hfc134ams*rdsigs
1073    
1074     pfcms = pfcms*rdsigs
1075    
1076     sf6ms = sf6ms*rdsigs
1077     ! ===
1078     #endif
1079    
1080     bcms = bcms*rdsigs
1081     ocms = ocms*rdsigs
1082    
1083     c 062295
1084     c h2o2ms = h2o2ms*rdsigs
1085     c
1086     !
1087     #endif
1088    
1089     BLDATA(I,J,4)=LMAX 6739.
1090     DO 8180 L=1,LMAX 6740.
1091     AJL(J,L,12)=AJL(J,L,12)+(THM-T(I,J,L))*PK(I,J,L)*P(I,J) 6741.
1092     T(I,J,L)=THM 6742.
1093    
1094     #if ( defined CPL_CHEM )
1095     !
1096     ! --- 03/23/95
1097     !
1098     cfc11(i,j,l) = cfc11ms
1099    
1100     cfc12(i,j,l) = cfc12ms
1101    
1102     xn2o(i,j,l) = xn2oms
1103    
1104     o3(i,j,l) = o3ms
1105    
1106     co(i,j,l) = coms
1107    
1108     zco2(i,j,l) = zco2ms
1109    
1110     xno(i,j,l) = xnoms
1111    
1112     xno2(i,j,l) = xno2ms
1113    
1114     xn2o5(i,j,l) = xn2o5ms
1115    
1116     hno3(i,j,l) = hno3ms
1117    
1118     ch4(i,j,l) = ch4ms
1119    
1120     ch2o(i,j,l) = ch2oms
1121    
1122     so2(i,j,l) = so2ms
1123    
1124     h2so4(i,j,l) = h2so4ms
1125    
1126     ! === if hfc, pfc, and sf6 are included:
1127     #ifdef INC_3GASES
1128     ! === 032698
1129     hfc134a(i,j,l) = hfc134ams
1130    
1131     pfc(i,j,l) = pfcms
1132    
1133     sf6(i,j,l) = sf6ms
1134     ! ===
1135     #endif
1136    
1137     bcarbon(i,j,l) = bcms
1138     ocarbon(i,j,l) = ocms
1139    
1140     c 062295
1141     c h2o2(i,j,l) = h2o2ms
1142     c
1143     !
1144     #endif
1145    
1146     8180 Q(I,J,L)=QMS 6743.
1147     IF(POLE) GO TO 8300 6744.
1148     C**** MIX MOMENTUM THROUGHOUT THE BOUNDARY LAYER AT NON-POLAR GRID BOXES6745.
1149     ID(1)=I+(J-1)*IM 6748.
1150     ID(2)=ID(1)+IM*JM*LM 6749.
1151     ID(3)=I+J*IM 6752.
1152     ID(4)=ID(3)+IM*JM*LM 6753.
1153     if(J.eq.JPR)then
1154     print *,'ID for J=',j
1155     print *,(ID(k),k=1,4)
1156     print *,'RA for J=',j
1157     print *,(RA(k),k=1,4)
1158     endif
1159     DO 8240 K=1,4 6754.
1160     UMS(K)=0. 6755.
1161     DO 8220 L=1,LMAX 6756.
1162     8220 UMS(K)=UMS(K)+UT(ID(K),1,L)*DSIG(L) 6757.
1163     8240 UMS(K)=UMS(K)*RDSIGS 6758.
1164     DO 8260 L=1,LMAX 6759.
1165     AJL(J,L,38)=AJL(J,L,38)+(UMS(1)-UT(I,J,L))*.5* 6760.
1166     * P(I,J)*RA(1) 6761.
1167     AJL(J+1,L,38)=AJL(J+1,L,38)+(UMS(3)- 6762.
1168     * UT(I,J+1,L))*P(I,J)*RA(3)*.5 6763.
1169     DO 8260 K=1,4 6764.
1170     if(J.eq.JPR)then
1171     print *,'L=',L,' K=',K
1172     print *,'ID(K)=',ID(K),' RA(K)=',RA(K)
1173     print *,'UMS(K)=',UMS(K),' UT(ID(K),1,L)=',UT(ID(K),1,L)
1174     endif
1175     8260 U(ID(K),1,L)=U(ID(K),1,L)+(UMS(K)-UT(ID(K),1,L))*RA(K) 6765.
1176     GO TO 8500 6766.
1177     C**** MIX MOMENTUM THROUGHOUT THE BOUNDARY LAYER AT POLAR GRID BOXES 6767.
1178     8300 JVPO=2 6768.
1179     IF(J.EQ.JM) JVPO=JM 6769.
1180     RAPO=2.*RAPVN(1) 6770.
1181     DO 8360 IPO=1,IM 6771.
1182     UMSPO=0. 6772.
1183     VMSPO=0. 6773.
1184     DO 8320 L=1,LMAX 6774.
1185     UMSPO=UMSPO+UT(IPO,JVPO,L)*DSIG(L) 6775.
1186     8320 VMSPO=VMSPO+VT(IPO,JVPO,L)*DSIG(L) 6776.
1187     UMSPO=UMSPO*RDSIGS 6777.
1188     VMSPO=VMSPO*RDSIGS 6778.
1189     DO 8340 L=1,LMAX 6779.
1190     U(IPO,JVPO,L)=U(IPO,JVPO,L)+(UMSPO-UT(IPO,JVPO,L))*RAPO 6780.
1191     V(IPO,JVPO,L)=V(IPO,JVPO,L)+(VMSPO-VT(IPO,JVPO,L))*RAPO 6781.
1192     8340 AJL(JVPO,L,38)=AJL(JVPO,L,38) 6782.
1193     * +(UMSPO-UT(IPO,JVPO,L))*P(1,J)*RAPO 6783.
1194     8360 CONTINUE 6784.
1195     8500 IM1=I 6793.
1196     do j=1,jm
1197     I=1
1198     if(J.eq.JPR.or.J.eq.-12)then
1199     print *,' J=',J,' after dry convection'
1200     print *,'U(I,J,1)=',U(I,J,1),' V(I,J,1)=',V(I,J,1)
1201     print *,'U(I,J+1,1)=',U(I,J+1,1),' V(I,J+1,1)=',V(I,J+1,1)
1202     endif
1203     enddo
1204     9000 CONTINUE 6794.
1205 jscott 1.3 ! print *,' From surf_ocean T2MZ'
1206     ! print *,T2MZ
1207 jscott 1.1 do 9001 J=1,JM
1208     c TSURFD(J)=TSURFD(J)+(BLDATA(1,J,2)-273.16)/24.
1209 jscott 1.3 TSURFD(J)=TSURFD(J)+(T2MZ(J)-273.16)/24.
1210 jscott 1.1 T2ML(J)=0.0
1211     9001 continue
1212     RETURN 6795.
1213     990 FORMAT ('0PPBL',3I4,14F8.2) 6818.
1214     991 FORMAT ('0SURFACE ',4I4,5F10.4,3F11.7) 6819.
1215     992 FORMAT ('0',I2,10F10.4/23X,4F10.4,10X,2F10.4/ 6820.
1216     * 33X,3F10.4,10X,2F10.4) 6821.
1217     993 FORMAT ('0',I2,10F10.4/23X,7F10.4/33X,7F10.4) 6822.
1218     994 FORMAT ('0',I2,11F10.4) 6823.
1219     END 6824.

  ViewVC Help
Powered by ViewVC 1.1.22