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

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

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


Revision 1.4 - (show annotations) (download)
Wed Sep 2 15:43:37 2009 UTC (15 years, 10 months ago) by jscott
Branch: MAIN
CVS Tags: HEAD
Changes since 1.3: +6 -1 lines
more flexible in allowing different time steps

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

  ViewVC Help
Powered by ViewVC 1.1.22