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

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

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


Revision 1.4 - (show annotations) (download)
Mon Apr 23 21:20:17 2007 UTC (18 years, 3 months ago) by jscott
Branch: MAIN
Changes since 1.3: +122 -71 lines
bring igsm atmos code up to date

1 c source sokolov users 76203 Apr 7 12:46 atmosphere.F
2
3 #include "ctrparam.h"
4
5 ! ==========================================================
6 !
7 ! Atmosphere.F: Former main program of the MIT Global Climate and
8 ! Biogeochemistry Model.
9 !
10 ! ----------------------------------------------------------
11 !
12 ! Repacking Note: This version is combined with main.f of
13 ! several originally separated model versions such as
14 ! CliChem 3.0, CliChemNem, as well as MODEL24x11.
15 !
16 ! The chemistry module, ocean CO2 module,
17 ! and TEM module are all controlled by cpp now.
18 !
19 ! Chien Wang
20 ! MIT Joint Program for the Science and Policy
21 ! of Global Change
22 !
23 ! ----------------------------------------------------------
24 !
25 ! Revision History:
26 !
27 ! When Who What
28 ! ---- ---------- -------
29 ! 073100 Chien Wang Created from CliChem3.0 and MODEL24x11
30 ! 100500 Andrei Converted to subroutine. When called first
31 ! time does only initialization.
32 ! Creates monthly data for TEM
33 !
34 ! ==========================================================
35
36
37 SUBROUTINE ATMOSPHERE(DTATM,mndriver)
38 C**** MD2G04 BD2G04 MD2G04 01/02/93 0.1
39 C**** OPT(3) 0.2
40 C**** 0.3
41 C**** Dynamics and physics programs for 2-D model. 0.4
42 C**** Like D2G04 but run on work station. 0.5
43
44 #include "BD2G04.COM" 1.
45 #include "ODIFF.COM"
46 #include "run.COM" 1.
47
48 #include "chem_para"
49 #include "chem_com"
50
51 #if ( defined CPL_CHEM )
52 !
53 #include "chem_tmp"
54 integer hrcnt, cnt3hr(8) ! for ozone impact
55 real sfc3hro3(nlat,8)
56 creal*4 sfc3hro3(nlat,8)
57 !
58 #endif
59
60 COMMON/INTA/COE1(01,01,01),COE2(01,01,01) 1.1
61 c DIMENSION ACO(36,01),BSI(36,01),CCO(36,01),DSI(36,01) 1.2
62 COMMON/SPEC1/ 1.3
63 * XA(IM0,JM0),XB(IM0,JM0),YA(IM0,JM0),YB(IM0,JM0),ZA(IM0,JM0)
64 &,ZB(IM0,JM0) 1.4
65 COMMON/SPEC2/KM,KINC,COEK,C3LAND(IO0,JM0),C3OICE(IO0,JM0) 1.5
66 * ,C3LICE(IO0,JM0),WMGE(IO0,JM0),TSSFC(1,JM0,4) 1.6
67 COMMON/EPARA/VTH(JM0,LM0),WTH(JM0,LM0),VU(JM0,LM0),VV(JM0,LM0)
68 &,DQSDT(JM0,LM0) 1.7
69 * ,DWV(JM0),PHIT(JM0,LM0),TPRIM2(JM0,LM0),WU(JM0,LM0),CKS,CKN 1.8
70 * ,WQ(JM0,LM0),VQ(JM0,LM0),MRCHT 1.9
71 COMMON U,V,T,P,Q 2.
72 COMMON/WORK1/WORKX(IM0,JM0,LM0),UT(IM0,JM0,LM0),VT(IM0,JM0,LM0), 3.
73 * TT(IM0,JM0,LM0),PT(IM0,JM0),QT(IM0,JM0,LM0) 4.
74 COMMON/WORK2/UX(IM0,JM0,LM0),VX(IM0,JM0,LM0),TX(IM0,JM0,LM0)
75 &,PX(IM0,JM0) 5.
76 COMMON/OLDZO/ZMLOLD(IO0,JM0)
77 C COMMON/KEYS/KEYNR(42,50) 8.
78 CHARACTER C*4,CYEAR*4,CMMND*80 8.1
79 DIMENSION C(39),JC(100),RC(161) 8.2
80 EQUIVALENCE (JC(1),IM),(C(1),XLABEL(1),LABEL1),(RC(1),TAU) 8.3
81 CHARACTER*8 LABSSW,LABEL1,OFFSSW/'XXXXXXXX'/ 8.4
82 LOGICAL EVENT,wr25,HPRNT,TRANSR,LHORDIF 9.
83 & ,CONTRR,OBSFOR,FIRST,NOCLM
84 common/conprn/HPRNT,JPR,LPR
85 common/wrcom/wr25,TRANSR,CONTRR,OBSFOR
86 COMMON/CWMG/WMGEA(JM0),NWMGEA(JM0),RIGA(JM0),DTAV(JM0),DQAV(JM0),
87 *Z0AV(JM0),WSAV(JM0),WS0AV(JM0),TAUAV(JM0)
88 COMMON/OCN/TG3M(1,JM0,12),RTGO(1,JM0,lmo),STG3(1,JM0),DTG3(1,JM0)
89 common/SURRAD/TRSURF(JM0,4),SRSURF(JM0,4)
90 dimension RTGOAV(JM0,lmo)
91 common/tprmtg/tprmg(JM0),ntprmg(JM0)
92 common/aexpc/AEXP,ISTRT1,ISTRTCHEM,LYEAREM
93 common/mixlr/Z1OAV(JM0),NZ1OAV(JM0)
94 common/flxio/FLIO(JM0),NFLIO(JM0)
95 common/surps/srps(JM0+3),nsrps
96 c character *19 buf
97 c character *8 buf1
98 character *120 file1,file2,plotfl,nwrfl,qffile,clfile
99 CHARACTER*4 AMONTH(12),JMONTHPR
100 common/files/file1,file2,flotfl,nwrfl,qffile,clfile
101 common/PRNT1/NCOMP
102 common/Dscale/DWAV0(JM0)
103 COMMON/CO2TRND/ALFFOR,CO2TR
104 COMMON/FRMIC/ FRMDICE(JM0)
105 common/ BACKGRGHG/GHGBGR(5)
106 DATA AMONTH/'JAN','FEB','MAR','APR','MAY','JUNE','JULY'
107 & ,'AUG','SEP','OCT','NOV','DEC'/
108 dimension fluxnep(jm0)
109 logical odifcarbon
110
111 #if ( defined CLM )
112 # include "CLM.h"
113 # if ( defined CPL_TEM )
114 C For TEM
115 # include "TEM.h"
116 # endif
117 #endif
118
119 #if ( defined CPL_OCEANCO2 )
120 #include "OCM.COM"
121 common /Garyflux/pC_atm(jm0),wind_amp,fluxco2(jm0)
122 # if ( defined ML_2D )
123 common/Garyclim/tggary(jm0),wsgary(jm0),areaml(jm0),arsrf(jm0)
124 common/Garydiff/depthml(jm0),edzcar(jm0),dzg(lmo),dzog(lmo-1),
125 &Rco2(jm0,lmo),edohd(lmo),zg(lmo),focean(jm0)
126 common /Garychem/Hg(jm0)
127 common/qfl/QFLUX(JM0,0:13),ZOAV(JM0),QFLUXT(JM0)
128 common /Garyvdif/iyearocm,vdfocm,acvdfc
129 common /Garyvlog/odifcarbon,ocarcont
130 common /Garykvct/cfkvct,edzcart(jm0)
131 # endif
132 #endif
133
134 INTEGER dtatm, mndriver !routine arguments jrs
135
136 #if ( defined OCEAN_3D || defined ML_2D )
137 #include "AGRID.h"
138 Cjrs elimated COM file/moved elsewhere#include "HRD4OCN.COM"
139 dimension oimeltt(jm0),dhdtav(jm0),devdtav(jm0)
140 #endif
141
142 C **** CLEAR SKY
143 common/clrsk/CLEAR(JM0),NCLR(JM0),AJCLR(JM0,12),BJCLR(JM0,12),
144 * CJCLR(JM0,12)
145 ! common/TSUR/TSURFC(JM0,0:13),TSURFT(JM0),TSURFD(JM0),DTSURF(JM0)
146 #include "TSRF.COM"
147 real TSURFW(JM0),TLANDW(JM0)
148 integer CLEAR
149
150 common /ATCO2/atm_co2(jm0)
151
152 #if ( defined CPL_NEM )
153 C For Emission
154 c === 031097
155 real ECH4COR(JM0),ECH4PY(JM0,12), ECH4OUT(JM0),EPJT(JM0)
156 &,ECH4CTR(JM0)
157 real EN2OCOR(JM0),EN2OPY(JM0,12), EN2OOUT(JM0),EPJTN2O(JM0)
158 &,EN2OCTR(JM0)
159 c ECH4CHIEN and EN2OCHIEN are used in chemistry model
160 common/EMFORCHIEN/ECH4CHIEN(JM0),EN2OCHIEN(JM0)
161 C For Emission
162 #endif
163 dimension NDAYMN(12)
164 data NDAYMN /31,28,31,30,31,30,31,31,30,31,30,31/
165 INTEGER JDOFM(13)
166 DATA JDOFM/0,31,59,90,120,151,181,212,243,273,304,334,365/
167 DATA FIRST/.TRUE./
168 SAVE
169
170 C AJCLR
171 C 1 SW INC AT P0 RD (AJ(1))
172 C 2 SW ABS BELOW P0 RD (AJ(2))
173 C 3 SW ABS BELOW P1 RD (AJ(3))
174 C 4 SW ABS AT Z0 RD (AJ(6))
175 C 5 SW INC AT Z0 RD (AJ(5))
176 C 6 LW INC AT Z0 RD (AJ(67))
177 C 7 NET LW AT Z0 SF (AJ(9))
178 C 8 NET LW AT P0 RD (AJ(7))
179 C 9 NET LW AT P1 RD (AJ(8))
180 C 10 NET RAD AT P0 DG (AJ(10))
181 C 11 NET RAD AT P1 DG (AJ(11))
182 C 12 NET RAD AT Z0 DG (AJ(12))
183 C **** CLEAR SKY
184 INTFX(XTAU)=INT(XTAU*XINT+.5) 10.
185 EVENT(XTAU)=MOD(ITAU,INTFX(XTAU)).LT.IDTHR 11.
186
187 !
188 ! --- assign input and output files
189 ! Note: Due to historical reasons, no all files are
190 ! assigned here - in case you want to search
191 ! something use
192 ! grep -i "needed characters" *.F
193 !
194 ! You have my sympathy.
195 !
196 ! Chien 080400
197 !
198
199 !#include "assign.inc"
200 c
201 if(FIRST) then
202 c
203 CALL CLOCKS (MNOW) 12.
204 C CALL ERRSET (206,1,0,0,0,301) 13.
205 C CALL ERRSET (208,256,-1) 14.
206 MBEGIN=MNOW 14.1
207 c CALL HARMO(36,1,24,ACO,BSI,CCO,DSI,1) 14.5
208 IPFLAG=0 15.
209 C CALL ENQJOB 16.
210 CALL INPUT 17.
211
212 #if ( defined CPL_CHEM )
213 !
214 ! --- Chemistry model
215 ! --- Set year and month index:
216 !
217 myyear = 1 !year from starting point
218 myyear = JYEAR-1976 !year from starting point
219 ! myyear = JYEAR-1891 !year from starting point
220 print *,'Emissioms for ',nchemyr,' year'
221 myyearlast = min(LYEAREM-1976,nchemyr) !last year of emission
222 ! myyearlast = min(LYEAREM-1891,nchemyr) !last year of emission
223 ! myyear = min(myyear,nchemyr)
224 myyear = min(myyear,myyearlast)
225 mymonth = 1 !month
226
227 ihaha = 1
228 ievenodd = 0 ! even hour 0
229 ! odd hour 1
230
231 call chembudget (p)
232 !
233 ! --- Set cfcnsf = 0.0
234 !
235 do k=1,nlev
236 cfcnsf(k) = 0.0
237 enddo
238 print *,'First year of emissions ', myyear
239 ! print *,'Emission will be fixed at year ',LYEAREM
240 print *,'Emission will be fixed at year ',1976+myyearlast
241 ! print *,'Emission will be fixed at year ',1891+myyearlast
242 !
243 #endif
244
245 C
246 KDISK0=500+KDISK
247 ndaa=3
248 c LHORDIF=.false.
249 LHORDIF=.true.
250 if(.not.LHORDIF)print *,' NO HOR. DIFFUSION for Q'
251 if(LHORDIF)print *,' HOR. DIFFUSION for Q after COND'
252 print *,' RADIATION EVERY ',NRAD/NDYN,' HOURS'
253
254 odifcarbon=.false.
255 #if ( defined CPL_OCEANCO2 && defined ML_2D )
256 odifcarbon=.true.
257 wind_amp=1.
258 dtco2=3600.*24.
259 c ncallgary=0
260 do j=1,jm
261 areaml(j)=dxyp(j)*(1-FDATA(1,J,2))
262 focean(j)=(1-FDATA(1,J,2))
263 DEPTHML(j)=ZOAV(j)
264 end do ! j
265 print *,' RCO2'
266 ! print 5001,((Rco2(j,k),j=1,jm),k=1,LMO)
267 print 5001,((Rco2(j,k)*1.e2,j=1,jm),k=1,LMO)
268 dzog(1)=10./SQRT(1.7010587)
269 dzg(2)=10.
270 do l=2,lmo-1
271 dzog(L)=dzog(L-1)*1.7010587
272 dzg(L+1)=dzg(L)*1.7010587
273 end do
274 zg(1)=50.
275 dzg(1)=100.
276 do l=2,lmo
277 zg(l)=zg(l-1)+0.5*(dzg(l-1)+dzg(l))
278 end do
279 do l=1,lmo
280 c edohd(l)=2.5e4/(zg(l)/zg(1))**1.0
281 C New coefficients for horizontal diffusion 11/16/00
282 edohd(l)=1.55e4-9.231e3*(atan((zg(l)-300)/50))
283 end do
284 #endif
285
286
287 #if ( defined CPL_CHEM )
288 !
289 ! --- Initialization of chemistry model:
290 !
291 c call cheminit(ISTRT1,T,q)
292 c 11/30/2000 ISTRTCHEM for restart of chemistry model
293 call cheminit(ISTRTCHEM,T,q)
294 c print *,'H2SO4 after cheminit ',h2so4(1,33,1)
295
296
297 !
298 ! --- tmp output
299 !
300 copen(124,file='OUTPUT/hro3',form='unformatted',
301 c & status='unknown')
302 hrcnt = 1
303 cnt3hr(:) = 0
304 sfc3hro3(:,:) = 0.0
305 !
306 #endif
307
308 print *,' IRAND=',IRAND
309 print *,' NCNDS=',NCNDS
310 print *,' after INPUT MRCHT=',MRCHT
311 JDAY00=JDAY-1
312 AEXP4=AEXP
313 TAU4=TAU
314 print *,' ISTRT1=',ISTRT1
315 if(ISTART.eq.2.or.ISTRT1.eq.0)then
316 nwr=1
317 WRITE (546) AEXP4,TAU4,XLABEL
318 write(547)AEXP,nwr
319 elseif(ISTART.eq.10)then
320 read(547)AEXPX,nwr
321 print *,' NWR=',nwr
322 if(abs(AEXPX-AEXP).gt.0.05)then
323 print *,' DISAGREEMENT BETWEEN AEXPX AND AEXP FILE 47'
324 print *,' AEXPX=',AEXPX,' AEXP=',AEXP
325 stop
326 endif
327 READ (546) AEXP4
328 if(abs(AEXP4-AEXP).gt.0.05)then
329 print *,' DISAGREEMENT BETWEEN AEXP4 AND AEXP FILE 46'
330 print *,' AEXP4=',AEXP4,' AEXP=',AEXP
331 stop
332 endif
333 C***
334 do 245 nr=1,nwr-1
335 READ (546)
336 245 continue
337 endif
338 WRITE(503) OFFSSW 17.1
339 REWIND 503 17.2
340 c CALL FRTR0(IO) 18.
341 KBGN=KINC+1 18.5
342 KM2=KM*2-1 18.51
343 KM3=KM2 18.52
344 KB1=3 18.53
345 KB2=5 18.54
346 KB3=5 18.55
347 IS=IM 18.56
348 IF(KM.EQ.1) IS=1 18.57
349 FIO=IO 18.58
350 JMM2=JM-2 18.59
351 HR24=24. 18.6
352 HR12=12. 18.61
353 MSTART=MNOW+MDYN+MCNDS+MRAD+MSURF+MDIAG+MELSE 19.
354 LMT3P1=LM*3+1 20.
355 C**** INITIALIZE TIME PARAMETERS 21.
356 DTHR=DT/3600. 22.
357 IDTHR=INTFX(DTHR) 23.
358 DTFS=DT*2./3. 24.
359 DTLF=2.*DT 25.
360 NDYNO=MOD(NDYN,2) 26.
361 I24=INTFX(24.) 27.
362 NSTEP0=.5+TAUI/DTHR 28.
363 NSTEP=INT(.5+TAU/DTHR)-NSTEP0 29.
364 NSTEP1=NSTEP 29.5
365 NSTEP2=NSTEP 29.6
366 MRCHT=0. 29.7
367 ITAU=(NSTEP+NSTEP0)*IDTHR 30.
368 TAU=FLOAT(ITAU)/XINT 31.
369 IDAY=1+ITAU/I24 32.
370 TOFDAY=(ITAU-(IDAY-1)*I24)/XINT 33.
371 ! if(ISTART.eq.2.or.ISTRT1.eq.0.and..not.CONTRR)then
372 ! do 458 j=1,JM
373 ! TSURFD(j)=0.
374 ! TSURFT(j)=0.
375 ! 458 continue
376 ! endif
377 if(JDATE.eq.100)then
378 print *,JDATE,JMONTH,JYEAR
379 print *,' main before daily0'
380 print *,' T1 ocean'
381 print 5001,(ODATA(1,j,1),j=1,JM0)
382 print *,' T2 ocean'
383 print 5001,(ODATA(1,j,4),j=1,JM0)
384 print *,' T3 ocean'
385 print 5001,(ODATA(1,j,5),j=1,JM0)
386 print *,' sea ice'
387 print 5002,(ODATA(1,j,2),j=1,JM0)
388 endif
389 CALL DAILY_NEW0 34.
390 print *,' Main after DAILYNEW0 JYEAR=',JYEAR
391 print *,"DT2MGL"
392 print *,DT2MGL
393 print *,"DT2MLD"
394 print *,DT2MLD
395 #if( !defined OCEAN_3D&& !defined ML_2D )
396 CALL DAILY_OCEAN
397 print *,' AFTER DAILY_OCEAN IDAY=',IDAY,' IYEAR=',IYEAR
398 print *,' JYEAR=',JYEAR,' JDAY=',JDAY
399 print *,' JDATE=',JDATE,' JMONTH=',JMONTH
400 #endif
401 if(JDATE.eq.100)then
402 print *,JDATE,JMONTH,JYEAR
403 print *,' main after daily0'
404 print *,' T1 ocean'
405 print 5001,(ODATA(1,j,1),j=1,JM0)
406 print *,' T2 ocean'
407 print 5001,(ODATA(1,j,4),j=1,JM0)
408 print *,' T3 ocean'
409 print 5001,(ODATA(1,j,5),j=1,JM0)
410 print *,' sea ice'
411 print 5002,(ODATA(1,j,2),j=1,JM0)
412 endif
413 c print *,' Z12O'
414 c print *,(Z12O(1,j),j=1,JM0)
415 c99 CONTINUE 34.993
416 CALL CLOCKS (MLAST) 35.
417 MINC=MNOW-MLAST 36.
418 MELSE=MELSE+MINC 37.
419 PERCNT=100.*MELSE/(MSTART-MLAST+1.E-5) 38.
420 c WRITE (6,901) IDAY,TOFDAY,JDATE,JMONTH,MINC,MELSE,PERCNT,TAU 39.
421 DOPK=1. 40.
422 MODD5K=1000 41.
423 IF(TAU.GE.TAUE) GO TO 820 42.
424 HPRNT=.TRUE.
425 HPRNT=.FALSE.
426 JPR=7
427 JPR=1
428 JPR=14
429 LPR=1
430 print *,' MAIN MRCHT=',MRCHT
431 c
432 c print *,' TAUE=',TAUE
433 TAUEM=TAUE
434 FIRST=.FALSE.
435 C
436 c TAU for coupler
437 c
438 TAUATM=TAU
439 MONTHATM=JMONTH
440 JDATEATM=JDATE
441 JYEARATM=JYEAR
442 C
443 #if ( defined CPL_OCEANCO2 )
444 do j=1,jm
445 fluxco2(j)=0.0
446 enddo
447 #endif
448 #if ( defined CPL_CHEM) && ( defined CPL_TEM )
449 C For TEM
450 if(ISTRT1.eq.0) then
451 c New run
452 c Reading from flin_nep
453 read(537)adupt,temco2
454 else
455 c Restart of the run
456 c Reading from last_nep
457 read(367)adupt,temco2
458 ! & ,temch4,temn2o
459 rewind 367
460 endif
461 !
462 ! adupt= 1.459814341652516
463 ! adupt= 0.9078891180588442
464 ! adupt= 0.25
465 ! adupt= -0.1123070421398009
466 !
467 ! adupt= adupt+0.9 ! for vs23
468
469 aduptd=adupt/(365.*JM)
470 temnepgl=0.0
471 do j=1,jm
472 temnepgl=temnepgl+temco2(j)
473 enddo
474 print *,'ADNEP=',adupt
475 print *,'Initial NEP=',adupt+temnepgl*1.e-3
476 temup0=0.0
477 #endif
478 #if (!defined OCEAN_3D && !defined ML_2D)
479 if(TRANSR)then
480 if(LMO.eq.11) then
481 call ODIFS
482 elseif(LMO.eq.12) then
483 call ODIFS12
484 else
485 Print *,' Wrong LMO',LMO
486 stop
487 endif
488 endif
489 #endif
490 !#if (defined PREDICTED_GASES)
491 #if (defined CPL_TEM || defined CPL_OCEANCO2 )
492 if(OBSFOR) then
493 call obsco2(iyear,imontha)
494 mnobco2=imonth
495 endif
496 #endif
497 !#endif
498 CJRS removed below from ocean_3d
499 #ifdef ML_2D
500 do j=1,jm
501 do i=1,io
502 CLAND4OCEAN(i,j)= C3LAND(I,J)
503 enddo
504 enddo
505 IDAYM=IDAY
506 JDAYM=JDAY
507 JDATEM=JDATE
508 JMONTHM=JMONTH
509 JYEARM=JYEAR
510 TAUML=TAU
511 TOFDAYML=TOFDAY
512 KOCEANM=KOCEAN
513 #endif
514 #if (defined OCEAN_3D || defined ML_2D)
515 do l=1,lm
516 sigfl(l)=sig(l)
517 enddo
518 print *,'SIGFL'
519 print *,sigfl
520 #endif
521 #if ( defined CPL_CHEM )
522 do j=1,jm
523 atm_co2(j)=zco2(1,j,1)
524 & *28.97296245/44.0*1.e-9
525 !ppb(m) to kg per volume base
526 & *1.e6
527 enddo
528 #else
529 if(.not.OBSFOR) then
530 CFF=1.
531 if (CO2.gt.0.0)CFF=CO2
532 do j=1,jm
533 atm_co2(j)=CFF*GHGBGR(1)
534 enddo
535 endif
536 #endif
537 #if (defined CPL_TEM || defined CPL_OCEANCO2 )
538 print *,'ATM_CO2'
539 print *,atm_co2
540 #endif
541 JDAYLAST=-1
542 ncallclm=0
543 NOCLM=.true.
544 #if ( defined CLM )
545 NOCLM=.false.
546 #endif
547 co2mmavd=0.0
548 print *,' atmosphere DTATM=',DTATM
549 print *,' It is running'
550 print *,'End of atmospheric model initialization'
551 print *,' '
552 print *,' '
553 print *,' JMONTHM= ',JMONTHM
554 print *,' TOFDAYML= ',TOFDAYML
555 print *,' '
556 return
557 endif ! first
558 C**** 43.
559 C**** MAIN LOOP 44.
560 C**** 45.
561 C
562 cprint *,' atmosphere TAU=',tau
563 TAUE=TAU+DTATM
564 c HPRNT=TAU.ge.17520.0.and.TAU.lt.17545.0
565 c print *,' TAUE=',TAUE
566 #if ( defined OCEAN_3D || defined ML_2D)
567 do j=1,jm0
568 tauu(j)=0.
569 tauv(j)=0.
570 precip(j)=0.
571 evao(j)=0.
572 evai(j)=0.
573 hfluxo(j)=0.
574 hfluxi(j)=0.
575 dhfodtg(j)=0.
576 devodtg(j)=0.
577 dhfidtg(j)=0.
578 devidtg(j)=0.
579 dhfodtgeq(j)=0.
580 devodtgeq(j)=0.
581 dhfidtgeq(j)=0.
582 devidtgeq(j)=0.
583 tempr(j)=0.
584 cjrs change var name to arunoff
585 arunoff(j)=0.
586 solarinc_ice(j)=0.
587 solarnet_ice(j)=0.
588 solarinc_ocean(j)=0.
589 solarnet_ocean(j)=0.
590 Cjrs why? surfpr(j)=0.
591 naveo(j)=0.
592 navei(j)=0.
593 navrad(j)=0.
594 navrado(j)=0.
595 c
596 ps4ocean(j)=0.
597 do l=1,lm
598 qyz4ocean(j,l)=0.
599 tyz4ocean(j,l)=0.
600 enddo
601 c
602 enddo
603 #endif
604 #if ( defined OCEAN_3D && defined CPL_OCEANCO2 )
605 ! SECDAY=24.*3600.
606 SECDAY=DTATM*3600.
607 co2mmav=0.0
608 do j=1,jmocean
609 co2mmav=co2mmav+mmco2flux(j)
610 enddo
611 print *,'CO2F form ocean TAU=',TAU,co2mmav*12.e-15
612 Cjrs *** this block needs attention?? ***
613 Cjrs fluxco2(1)=fluxco2(1)+SECDAY*mmco2flux(1)
614 fluxco2(1)=fluxco2(1) + SECDAY*mmco2flux(2)
615 do j=2,jm-1
616 fluxco2(j)=fluxco2(j)+SECDAY*mmco2flux(j-1)
617 enddo
618 Cjrs fluxco2(JM)=fluxco2(JM)+SECDAY*mmco2flux(JMOCEAN)
619 fluxco2(JM)=fluxco2(JM) + SECDAY*mmco2flux(JM-1)
620 # endif
621
622 WLMMAX=0.0
623 C
624 100 IF(.NOT.EVENT(TAUT)) GO TO 200 46.
625 c HPRNT=TAU.ge.17520.00
626 NSTEP1=NSTEP 46.5
627 C**** WRITE RESTART INFORMATION ONTO DISK 47.
628 120 CALL RFINAL (IRAND) 48.
629 IF(NSTEP.EQ.NSTEP2) GO TO 116 48.3
630 DO 115 K=1,22 48.5
631 DO 115 J=1,JM 48.6
632 DO 115 I=2,IO 48.7
633 115 GDATA(I,J,K)=GDATA(1,J,K) 48.9
634 116 CONTINUE 48.91
635 if(wr25.and.ISTART.eq.2)then
636 print *,' main write'
637 print *,' T1 ocean'
638 print 5001,(ODATA(1,j,1),j=1,JM0)
639 print *,' T2 ocean'
640 print 5001,(ODATA(1,j,4),j=1,JM0)
641 print *,' T3 ocean'
642 print 5001,(ODATA(1,j,5),j=1,JM0)
643 REWIND KDISK0 49.
644 if(TRANSR)then
645 WRITE(KDISK0) AEXP,TAU,JC,C,RC,KEYNR,U,V,T,P,Q,ODATA,GDATA,BLDATA, 50.
646 * RQT,SRHR,TRHR,(AJ(K,1),K=1,KACC),TAU,TSSFC,CKS,CKN,WMGE,TPRIM2 51.
647 * ,MRCHT,TRSURF,SRSURF,TLANDW,TSURFW,DWAV0
648 * ,TG3M,RTGO,STG3,DTG3
649 print *,' STG3'
650 print 5001,(STG3(1,j),j=1,JM0)
651 print *,' DTG3/356'
652 print 5001,(DTG3(1,j)/365.,j=1,JM0)
653 print *,' RTGO'
654 print 5001,((RTGO(1,j,k),j=1,JM0),k=1,lmo)
655 5001 format(23f6.1)
656 5002 format(23f6.3)
657 else
658 WRITE(KDISK0) AEXP,TAU,JC,C,RC,KEYNR,U,V,T,P,Q,ODATA,GDATA,BLDATA, 50.
659 * RQT,SRHR,TRHR,(AJ(K,1),K=1,KACC),TAU,TSSFC,CKS,CKN,WMGE,TPRIM2 51.
660 * ,MRCHT,TRSURF,SRSURF,TLANDW,TSURFW,DWAV0
661 endif
662 REWIND KDISK0 52.
663 end if ! ISTART.eq.2
664 CALL CLOCKS (MNOW) 53.
665 MINC=MLAST-MNOW 54.
666 MELSE=MELSE+MINC 55.
667 PERCNT=100.*MELSE/(MSTART-MNOW+1.E-5) 56.
668 MLAST=MNOW 59.
669 C**** TEST FOR TERMINATION OF RUN 60.
670 200 READ (503,END=210) LABSSW 61.
671 c HPRNT=TAU.gt.45.0.and.TAU.lt.60.0
672 c HPRNT=TAU.gt.470.0.and.TAU.lt.550.0
673 NCOMP=0
674 210 REWIND 503 61.1
675 IF(LABSSW.EQ.LABEL1) KSS6=1 61.2
676 IF(KSS6.EQ.1) GO TO 800 62.
677 IF(TAU+.06125.GE.TAUE) GO TO 820 63.
678 JDAY00=JDAY
679 C**** IF TIME TO ZERO OUT DIAGNOSTIC ACCUMULATING ARRAYS, DO SO 64.
680 C**** (NORMALLY DONE AT THE BEGINNING OF A MONTH) 65.
681 IF (TAU.EQ.TAUI) GO TO 260 66.
682 IF(.NOT.EVENT(24.)) GO TO 290 67.
683 DO 250 K=1,13 68.
684 IF(JDAY.EQ.NDZERO(K)) GO TO 260 69.
685 250 CONTINUE 70.
686 GO TO 290 71.
687 260 CONTINUE
688 TAU0=TAU 72.
689 IDAY0=IDAY 73.
690 TOFDY0=TOFDAY 74.
691 JDATE0=JDATE 75.
692 JMNTH0=JMONTH 76.
693 JYEAR0=JYEAR 77.
694 DO 270 I=1,12 78.
695 270 IDACC(I)=0 79.
696 NODIFS=0
697 nsrps=0.
698 DO 280 K=1,KACC 80.
699 280 AJ(K,1)=0. 81.
700 do 5280 j=1,JM
701 tprmg(j)=0.
702 ntprmg(j)=0.
703 Z1OAV(j)=0.
704 NZ1OAV(j)=0
705 NFLIO(J)=0
706 FLIO(J)=0.
707 NCLR(J)=0
708 do 5282 k=1,lmo
709 RTGOAV(j,k)=0.
710 5282 continue
711 do 5281 n=1,12
712 AJCLR(J,n)=0.
713 BJCLR(J,n)=0.
714 CJCLR(J,n)=0.
715 5281 continue
716 5280 continue
717 do 5286 j=1,jm+3
718 srps(j)=0.
719 5286 continue
720
721 c#if ( defined CPL_OCEANCO2 && defined ML_2D )
722 c call zerogary(ncallgary)
723 c call zerogary
724 c#endif
725
726 CALL DIAG9A (1) 82.
727
728 #if ( defined CPL_NEM )
729 C For Emission
730 c === 031097
731 DO 5655 MONTH=1,12
732 IF(JDAY.LE.JDOFM(MONTH+1)) GO TO 5656
733 5655 CONTINUE
734 5656 MNHTEM=MONTH-1
735 if(MNHTEM.eq.0)MNHTEM=12
736 do J=1,JM
737 c ECH4CHIEN(J)=temch4(J)/NDAYMN(MNHTEM)/1000.
738 c EN2OCHIEN(J)=temn2o(J)/NDAYMN(MNHTEM)/1000.
739
740 ECH4CHIEN(J)=temch4(J)/NDAYMN(MNHTEM)
741 EN2OCHIEN(J)=temn2o(J)/NDAYMN(MNHTEM)
742
743 enddo
744
745 C For Emission
746 #endif
747
748 290 CONTINUE 83.
749 C**** 84.
750 C**** INTEGRATE DYNAMIC TERMS 85.
751 if(HPRNT)then
752 print *,' main before comp1 1',' TAU=',TAU,JDATE,JMONTH
753 #include "PRNT.COM"
754 print *,' GDATA(1,7,5)=',GDATA(1,7,5),' GDATA(1,7,6)='
755 * ,GDATA(1,7,6)
756 endif
757 C**** 86.
758 MODD5D=MOD(NSTEP,NDA5D) 87.
759 IF(MODD5D.EQ.0) CALL DIAG5A (2,0) 88.
760 c IF(NDYNO.EQ.1) then
761 c print *,NDYNO
762 c print *,jm,im,LMT3P1
763 c endif
764 DO 310 J=1,JM 89.
765 DO 300 L=1,LMT3P1 90.
766 DO 300 I=1,IM 91.
767 UX(I,J,L)=U(I,J,L) 92.
768 300 UT(I,J,L)=U(I,J,L) 93.
769 DO 310 L=1,LM 94.
770 DO 310 I=1,IM 95.
771 310 QT(I,J,L)=Q(I,J,L) 96.
772 if(HPRNT)then
773 print *,' main before comp1 2',' TAU=',TAU
774 #include "PRNT.COM"
775 print *,' GDATA(1,7,5)=',GDATA(1,7,5),' GDATA(1,7,6)='
776 * ,GDATA(1,7,6)
777 endif
778 C**** INITIAL FORWARD STEP, QX = Q + .667*DT*F(Q) 97.
779 NS=0 98.
780 MRCH=0 99.
781 CALL COMP1 (UX,VX,TX,PX,Q,U,V,T,P,Q,DTFS,NS) 100.
782 if(HPRNT)then
783 print *,' main after comp1',' TAU=',TAU,' MRCH=',MRCH
784 #include "PRNT.COM"
785 endif
786 C
787 IF(NDYNO.EQ.1) GO TO 320 101.
788 C**** INITIAL BACKWARD STEP IS ODD, QT = QT + DT*F(QX) 102.
789 MRCH=-1 103.
790 CALL COMP1 (UT,VT,TT,PT,QT,UX,VX,TX,PX,Q,DT,NS) 104.
791 if(HPRNT)then
792 print *,' main after comp1',' TAU=',TAU,' MRCH=',MRCH
793 #include "PRNT.COM"
794 endif
795 C
796 GO TO 360 105.
797 C**** INITIAL BACKWARD STEP IS EVEN, Q = Q + DT*F(QX) 106.
798 320 NS=1 107.
799 MODD5K=MOD(NSTEP+NS-NDYN+NDA5K,NDA5K) 108.
800 MRCH=1 109.
801 CALL COMP1 (U,V,T,P,Q,UX,VX,TX,PX,QT,DT,NS) 110.
802 if(HPRNT)then
803 print *,' main after comp1',' TAU=',TAU,' MRCH=',MRCH
804 #include "PRNT.COM"
805 endif
806 C
807 CD DIAGA SHOULD BE CALLED HERE BUT THEN ARRAYS MUST BE CHANGED 111.
808 c
809 c
810 C**** ODD LEAP FROG STEP, QT = QT + 2*DT*F(Q) 112.
811 340 MRCH=-2 113.
812 CALL COMP1 (UT,VT,TT,PT,QT,U,V,T,P,Q,DTLF,NS) 114.
813 if(HPRNT)then
814 print *,' main after comp1',' TAU=',TAU,' MRCH=',MRCH
815 #include "PRNT.COM"
816 endif
817 C
818 C**** EVEN LEAP FROG STEP, Q = Q + 2*DT*F(QT) 115.
819 360 NS=NS+2 116.
820 MODD5K=MOD(NSTEP+NS-NDYN+NDA5K,NDA5K) 117.
821 MRCH=2 118.
822 CALL COMP1 (U,V,T,P,Q,UT,VT,TT,PT,QT,DTLF,NS) 119.
823 if(HPRNT)then
824 print *,' main after comp1',' TAU=',TAU,' MRCH=',MRCH
825 #include "PRNT.COM"
826 endif
827 C
828 IF(NS.LT.NDYN) GO TO 340 122.
829 c IF(MOD(NSTEP+NS-NDYN+NDAA,NDAA).LT.MRCH) THEN 120.
830 CALL DIAGA (UT,VT,TT,PT,QT,NOCLM) 121.
831 if(HPRNT)then
832 print *,' main after DIAGA',' TAU=',TAU,' MRCH=',MRCH
833 #include "PRNT.COM"
834 endif
835 c ENDIF
836 c IF(NS.LT.NDYN) GO TO 340 122.
837 DOPK=1. 123.
838 CALL CLOCKS (MNOW) 124.
839 MINC=MLAST-MNOW 125.
840 MDYN=MDYN+MINC 126.
841 MLAST=MNOW 127.
842 PERCNT=100.*MDYN/(MSTART-MNOW+1.E-5) 128.
843 C 130.
844 CALL DIAG9A (2) 131.
845 C**** 133.
846 C**** INTEGRATE SOURCE TERMS 134.
847 C**** 135.
848 #if (!defined PREDICTED_GASES)
849 #if (defined CPL_TEM || defined CPL_OCEANCO2 )
850 if(OBSFOR) then
851 if(JMONTH.ne.AMONTH(mnobco2)) then
852 mnobco2=mnobco2+1
853 if(mnobco2.eq.13)mnobco2=1
854 call obsco2(jyear,mnobco2)
855 endif
856 endif
857 #endif
858 #endif
859 C
860 C
861 MODRD=MOD(NSTEP,NRAD) 136.
862 MODD5S=MOD(NSTEP,NDA5S) 137.
863 IF(MODD5S.EQ.0) IDACC(8)=IDACC(8)+1 138.
864 C 139.
865 C**** CONDENSTATION, SUPER SATURATION AND MOIST CONVECTION 140.
866 if(HPRNT)then
867 print *,' main before conds',' TAU=',TAU
868 #include "PRNT.COM"
869 endif
870
871
872 #if ( defined CPL_CHEM )
873 !
874 ! ===== Calculate airmass at grid point
875 ! ===== Chien Wang, 092395
876 !
877
878 call chemairmass(p4chem0)
879
880 !
881 ! === Calculating total mass of tracers
882 ! === with long residence times:
883 !
884 call chemmass1(cfc11, cfc11mass)
885 call chemmass1(cfc12, cfc12mass)
886 call chemmass1(xn2o, xn2omass)
887 call chemmass1(zco2, zco2mass)
888 call chemmass1(ch4, ch4mass)
889
890 ! === if hfc, pfc, and sf6 are included:
891 #ifdef INC_3GASES
892 ! === 032698
893 call chemmass1(hfc134a, hfc134amass)
894 call chemmass1(pfc, pfcmass)
895 call chemmass1(sf6, sf6mass)
896 ! ===
897 #endif
898
899 !
900 ! === Calculating advection and eddy diffusion:
901 !
902 call chemadv0 (dt)
903
904 !
905 ! === Calculate total n-s transport amount
906 ! === of cfc11 - temperary:
907 !
908 dth = 3600.0 ! for relexible setting of dt dt * 3.0
909 ! call chemtmp1 (dth,airmass0,p,pvv,cfc11)
910
911 ! === Readjust mass of tracers, 1:
912 !
913 ! -----------------------------
914 ! Use tropospheric life time (yr) to calculate mass
915 ! loss of tracers and at the same time compensate
916 ! all the numerical loss back to tracer's mass
917 ! which equavalent to use adjcoe = 1.0 for
918 ! chemmass2.f
919 !
920 ! Chien Wang, September 12,1995
921 ! ------------------------------
922
923 ! === 092595 update p
924
925 call chemairmass(p4chem1)
926
927 call chemmass6(46.0, 1.0, cfc11,cfc11mass)
928 call chemmass6(120.0,1.0, cfc12,cfc12mass)
929
930 ! === 102596
931 ! === close tau type of ocean uptake of co2:
932
933 call chemmass66(1.0, 1.0,zco2,zco2mass)
934
935 call chemmass6(150.0,1.0,xn2o,xn2omass)
936 call chemmass2(1.0,ch4, ch4mass )
937
938 ! === if hfc, pfc, and sf6 are included:
939 #ifdef INC_3GASES
940 ! === 032698
941 ! call chemmass6(14.6, 1.0,hfc134a, hfc134amass)
942 call chemmass2(1.0, hfc134a, hfc134amass)
943 call chemmass6(10000.0,1.0,pfc,pfcmass)
944 call chemmass6(3200.0, 1.0,sf6,sf6mass)
945 ! ===
946 #endif
947
948 !
949 ! === Calculate tropospheric gaseous reactions
950 ! === every nhr_for_chem hours:
951 !
952 dt_chem_h = dth*float(nhr_for_chem)
953 if(ievenodd.eq.0) then
954 call chemtrop0(0, T, q, dt_chem_h, 1)
955
956 c print *,'H2SO4 after chemtrop0 ',h2so4(1,33,1)
957 c print *,'SVIOD after chemtrop0 ',SVIOD(1,33,1)
958
959 !
960 ! --- tmp output
961 !
962 cnt3hr(hrcnt) = cnt3hr(hrcnt) + 1
963 sfc3hro3(1:nlat,hrcnt) = sfc3hro3(1:nlat,hrcnt)
964 & + (o3(1,1:nlat,1))*29.0/48.0
965
966 hrcnt = hrcnt + 1
967 if (hrcnt .gt. 8 ) hrcnt = 1
968 end if
969
970 ievenodd = ievenodd + 1
971 if(ievenodd.ge.nhr_for_chem) ievenodd = 0
972
973 !
974 ! === Calculating stratospheric processes:
975 ! === 092595
976 ! === adjust startospheric loss to whole global:
977 !
978 call chemmass1(cfc11, cfc11mass)
979 call chemmass1(cfc12, cfc12mass)
980 call chemmass1(xn2o, xn2omass)
981
982 call chemstrat (dt)
983
984 call chemmass2(1.0,cfc11,cfc11mass)
985 call chemmass2(1.0,cfc12,cfc12mass)
986 call chemmass2(1.0,xn2o, xn2omass)
987
988 !
989 ! === Get total mass of chemically active species:
990 !
991 ! call chemairmass(p)
992 ! write(6,*)"P before 2nd", p
993 c print *,'H2SO4 after chemairmass ',h2so4(1,33,1)
994
995 call chemmass1(cfc11, cfc11mass)
996 call chemmass1(cfc12, cfc12mass)
997 call chemmass1(xn2o, xn2omass)
998 call chemmass1(zco2, zco2mass)
999 call chemmass1(ch4, ch4mass)
1000
1001 ! === if hfc, pfc, and sf6 are included:
1002 #ifdef INC_3GASES
1003 ! === 032698
1004 call chemmass1(hfc134a, hfc134amass)
1005 call chemmass1(pfc, pfcmass)
1006 call chemmass1(sf6, sf6mass)
1007 ! ===
1008 #endif
1009
1010 #endif
1011
1012 c print *,'H2SO4 before CONDSE ',h2so4(1,33,1)
1013 CALL CONDSE(mndriver) 141.
1014 c print *,'H2SO4 after CONDSE ',h2so4(1,33,1)
1015
1016 #if ( defined CPL_CHEM )
1017 !
1018 ! === Calculate emission once per hour:
1019 !
1020 ! timeinhr=1./(365.*24.) !hourly emission
1021 ! call chememission(timeinhr)
1022 !
1023 !
1024 ! === Print hourly
1025 !
1026 ! call chemprt
1027 !
1028 ! ==============================================
1029 #endif
1030
1031 if(HPRNT)then
1032 print *,' main after conds',' TAU=',TAU
1033 #include "PRNT.COM"
1034 endif
1035 #if ( !defined CLM )
1036 CALL PRECIP_LAND(mndriver) 142.
1037 if(HPRNT)then
1038 print *,' main after preci',' TAU=',TAU
1039 #include "PRNT.COM"
1040 endif
1041 #endif
1042 CALL CLOCKS (MNOW) 143.
1043 MINC=MLAST-MNOW 144.
1044 MCNDS=MCNDS+MLAST-MNOW 145.
1045 MLAST=MNOW 146.
1046 C 147.
1047 CALL DIAG9A (3) 148.
1048 C**** RADIATION, SOLAR AND THERMAL 149.
1049 if(HPRNT)then
1050 print *,' main before radia',' TAU=',TAU
1051 #include "PRNT.COM"
1052 print *,' IRAND=',IRAND
1053 endif
1054 #if ( defined PREDICTED_GASES || defined PREDICTED_AEROSOL )
1055 call radia_chem
1056 #else
1057 if(OBSFOR) then
1058 CALL RADIAGSO
1059 else
1060 CALL RADIA
1061 endif
1062 #endif
1063 if(HPRNT)then
1064 print *,' main after radia',' TAU=',TAU
1065 #include "PRNT.COM"
1066 endif
1067 CALL CLOCKS (MNOW) 151.
1068 MINC=MINC+MLAST-MNOW 152.
1069 MRAD=MRAD+MLAST-MNOW 153.
1070 MLAST=MNOW 154.
1071 if(HPRNT)then
1072 print *,' main before diag9a',' TAU=',TAU
1073 #include "PRNT.COM"
1074 endif
1075 C 155.
1076 CALL DIAG9A (4) 156.
1077 if(HPRNT)then
1078 print *,' main after diag9a',' TAU=',TAU
1079 #include "PRNT.COM"
1080 endif
1081 C**** SURFACE INTERACTION AND GROUND CALCULATION 157.
1082
1083 #if ( defined CLM )
1084 if(HPRNT)then
1085 print *,' main before surf4clm',' TAU=',TAU
1086 #include "PRNT.COM"
1087 endif
1088 CALL SUR4CLM
1089 if(HPRNT)then
1090 print *,' main after surf4clm',' TAU=',TAU
1091 #include "PRNT.COM"
1092 endif
1093 i=1
1094 do j=1,jm
1095 pcpl4clm(i,j)=pcpl4clm(i,j)*prlnd2total(j,mndriver)
1096 pcpc4clm(i,j)=pcpc4clm(i,j)*prlnd2total(j,mndriver)
1097 enddo
1098 ! print *,' main after surf4clm',' TAU=',TAU
1099 ! print ('2(12f7.2,/,11f7.2,/)'),ps4clm,pcpl4clm,
1100 ! & pcpc4clm,tpr4clm,
1101 ! & tsl4clm,
1102 ! & qs4clm,ws4clm
1103 ! & ,us4clm,vs4clm,
1104 ! & dsw4clm,
1105 ! & dlw4clm,pco24clm
1106 ! & ,swinr4clm,swvis4clm
1107 #if ( defined DATA4TEM )
1108 c if(JYEAR.gt.20)then
1109 write (935),ps4clm,pcpl4clm,
1110 & pcpc4clm,tpr4clm,
1111 & tsl4clm,
1112 & qs4clm,ws4clm
1113 & ,us4clm,vs4clm,
1114 & dsw4clm,
1115 & dlw4clm,pco24clm
1116 & ,swinr4clm,swvis4clm
1117 c endif
1118 #endif
1119
1120 ncallclm=ncallclm+1
1121 c print *,'before clm4mit2d ncallclm=',ncallclm
1122 if(HPRNT)then
1123 print *,' main before clm4mit2d',' TAU=',TAU
1124 #include "PRNT.COM"
1125 endif
1126 call clm4mit2d
1127 if(HPRNT)then
1128 print *,' main after clm4mit2d',' TAU=',TAU
1129 #include "PRNT.COM"
1130 endif
1131
1132 c if(JYEAR.gt.20)then
1133 ! write (934),tau,tsoiclm,snwdclm,snwcclm,
1134 ! & lwuclm,tref2mclm,tflxclm,tgndclm,
1135 ! & lhfclm,shfclm,tauxclm,tauyclm,
1136 ! & asdirclm,aldirclm,asdifclm,aldifclm,
1137 ! & sroclm,ssrclm,glrclm
1138 ! &,h2olclm,h2oiclm
1139 c endif
1140 ! print *,' main after clm4mit2d',' TAU=',TAU
1141 ! print ('2(12f7.2,/,11f7.2,/)'),tsoiclm,snwdclm,snwcclm,
1142 ! & lwuclm,tref2mclm,tflxclm,tgndclm,
1143 ! & lhfclm,shfclm,tauxclm,tauyclm,
1144 ! & asdirclm,aldirclm,asdifclm,aldifclm
1145
1146 CALL SURF_CLM
1147 CALL SURF_OCEAN
1148 CALL GR_CLM
1149 if(HPRNT)then
1150 print *,' main after surfc',' TAU=',TAU
1151 #include "PRNT.COM"
1152 endif
1153 #else
1154 ! CALL SURF_LAND
1155 ! CALL SURF_OCEAN
1156 CALL SURFCE ! 07/14/2006
1157 if(HPRNT)then
1158 print *,' main after surfc',' TAU=',TAU
1159 #include "PRNT.COM"
1160 endif
1161 CALL GRLAND 159.
1162 #endif
1163
1164 #if ( defined OCEAN_3D || defined ML_2D)
1165 CALL GRFOROCEAN
1166 #else
1167 CALL GROCEAN(mndriver) 159.
1168 #endif
1169 c
1170 if(HPRNT)then
1171 print *,' main after groun',' TAU=',TAU
1172 #include "PRNT.COM"
1173 endif
1174 c print *,'H2SO4 before DRYCNV ',h2so4(1,33,1)
1175 CALL DRYCNV 160.
1176 c print *,'H2SO4 after DRYCNV ',h2so4(1,33,1)
1177 if(HPRNT)then
1178 print *,' main after drycn',' TAU=',TAU
1179 #include "PRNT.COM"
1180 endif
1181 CALL CLOCKS (MNOW) 161.
1182 MINC=MINC+MLAST-MNOW 162.
1183 MSURF=MSURF+MLAST-MNOW 163.
1184 MLAST=MNOW 164.
1185 CALL DIAG9A (5) 165.
1186 C**** STRATOSPHERIC MOMENTUM DRAG 166.
1187 CALL SDRAG(WLMMAX,JWLMMAX) 167.
1188 if(HPRNT)then
1189 print *,' main after sdrag',' TAU=',TAU
1190 #include "PRNT.COM"
1191 endif
1192 CALL CLOCKS (MNOW) 168.
1193 MINC=MINC+MLAST-MNOW 169.
1194 MSURF=MSURF+MLAST-MNOW 170.
1195 MLAST=MNOW 171.
1196 C 172.
1197 CALL DIAG9A (6) 173.
1198 MSRCE=MCNDS+MRAD+MSURF 174.
1199 PERCNT=100.*MSRCE/(MSTART-MNOW+1.E-5) 175.
1200 C**** SEA LEVEL PRESSURE FILTER 177.
1201 NFILTR= NDYN 177.5
1202 IF(MFILTR.LE.0.OR.MOD(NSTEP,NFILTR).NE.0) GO TO 500 178.
1203 IDACC(10)=IDACC(10)+1 179.
1204 if(HPRNT)then
1205 print *,' main before filtr',' TAU=',TAU
1206 #include "PRNT.COM"
1207 endif
1208 C
1209 C 180.
1210 C ****************
1211 if(LHORDIF)then
1212 DTDIF=3600.
1213 if(JM.eq.24)then
1214 CALL HORDIFF(DTDIF)
1215 else if(JM.eq.46)then
1216 CALL HORDIFFALL(DTDIF)
1217 else
1218 print *,' Wromg JM=',JM
1219 stop
1220 endif
1221 end if
1222 C ****************
1223 c CALL FILTER 181.
1224 C
1225 #if ( defined CPL_CHEM )
1226 !
1227 ! === Readjust total mass of tracers:
1228 !
1229 call chemairmass(p)
1230
1231 call chemmass2(1.00,cfc11,cfc11mass)
1232 call chemmass2(1.00,cfc12,cfc12mass)
1233 call chemmass2(1.00,xn2o ,xn2omass )
1234 call chemmass2(1.00,zco2,zco2mass)
1235 call chemmass2(1.00,ch4, ch4mass )
1236
1237 ! === if hfc, pfc, and sf6 are included:
1238 #ifdef INC_3GASES
1239 ! === 032698
1240 call chemmass2(1.0,hfc134a, hfc134amass)
1241 call chemmass2(1.0,pfc, pfcmass)
1242 call chemmass2(1.0,sf6, sf6mass)
1243 ! ===
1244 #endif
1245
1246 !
1247 ! === Accumulative calculation prepared for
1248 ! === carrying out monthly average:
1249 !
1250 call chemmonth1
1251 !
1252 #endif
1253
1254 CALL CLOCKS (MNOW) 182.
1255 MDYN=MDYN+MLAST-MNOW 183.
1256 MLAST=MNOW 184.
1257 C 185.
1258 CALL DIAG9A (7) 186.
1259 C**** 187.
1260 C**** UPDATE MODEL TIME AND CALL DAILY IF REQUIRED 188.
1261 C**** 189.
1262 500 NSTEP=NSTEP+NDYN 190.
1263 ITAU=(NSTEP+NSTEP0)*IDTHR 191.
1264 TAU=FLOAT(ITAU)/XINT 192.
1265 IDAY=1+ITAU/I24 193.
1266 TOFDAYPR=TOFDAY+1.00
1267 TOFDAY=(ITAU-(IDAY-1)*I24)/XINT 194.
1268 IF(.NOT.EVENT(24.)) GO TO 550 195.
1269 C 196.
1270 do J=1,JM0
1271 TSURFW(J)=TSURFD(J)
1272 TLANDW(J)=TLANDD(J)
1273 enddo
1274
1275 JDATECLM=JDATE
1276 JDATEPR=JDATE
1277 JMONTHPR=JMONTH
1278 JYEARPR=JYEAR
1279 CALL DAILY_NEW 197.
1280 c print *,' AFTER DAILY_NEW IDAY=',IDAY,' IYEAR=',IYEAR
1281 c print *,' JYEAR=',JYEAR,' JDAY=',JDAY
1282 #if( !defined OCEAN_3D && !defined ML_2D )
1283 CALL DAILY_OCEAN
1284 c print *,' AFTER DAILY_OCEAN IDAY=',IDAY,' IYEAR=',IYEAR
1285 c print *,' JYEAR=',JYEAR,' JDAY=',JDAY
1286 c print *,' JDATE=',JDATE,' JMONTH=',JMONTH
1287 #endif
1288 if(JDATE.eq.100)then
1289 print *,JDATE,JMONTH,JYEAR
1290 print *,' main after daily'
1291 print *,' T1 ocean'
1292 print 5001,(ODATA(1,j,1),j=1,JM0)
1293 print *,' T2 ocean'
1294 print 5001,(ODATA(1,j,4),j=1,JM0)
1295 print *,' T3 ocean'
1296 print 5001,(ODATA(1,j,5),j=1,JM0)
1297 print *,' sea ice'
1298 print 5002,(ODATA(1,j,2),j=1,JM0)
1299 endif
1300 CALL CLOCKS (MNOW) 198.
1301 MELSE=MELSE+(MLAST-MNOW) 199.
1302 MLAST=MNOW 200.
1303 NDAILY=SDAY/DT 201.
1304
1305 #if ( defined CPL_CHEM )
1306 !
1307 ! === Calculate air mass, second step:
1308 !
1309 i=1
1310 call chemairmass(p)
1311
1312 !
1313 ! === Calculate emission once per day:
1314 !
1315 c zxy=0.
1316 c do j=1,jm
1317 c zxy=zxy+zco2(1,j,1)
1318 c & *28.97296245/44.0*1.e-3
1319 c enddo
1320 c print *,' CO2 before emission ',zxy/jm
1321
1322 timeinday=1./(365.) !daily emission
1323 call chememission(timeinday)
1324
1325 c zxy=0.
1326 c do j=1,jm
1327 c zxy=zxy+zco2(1,j,1)
1328 c & *28.97296245/44.0*1.e-3
1329 c enddo
1330 c print *,' CO2 after emission ',zxy/jm
1331 !
1332 #endif
1333
1334 C 202.
1335 CALL DIAG9A (8) 203.
1336 #if ( !defined OCEAN_3D && !defined ML_2D )
1337 IF(KOCEAN.EQ.1) THEN
1338 DO 540 J=1,JM 203.11
1339 DO 540 I=1,IM 203.12
1340 AIJ(I,J,59)=AIJ(I,J,59)+ODATA(I,J,4) 203.13
1341 540 AIJ(I,J,60)=AIJ(I,J,60)+ODATA(I,J,5) 203.14
1342 if(TRANSR)then
1343 if(LMO.eq.11) then
1344 CALL ODIFS
1345 elseif(LMO.eq.12) then
1346 CALL ODIFS12
1347 else
1348 Print *,' Wromng LMO',LMO
1349 stop
1350 endif
1351 NODIFS=NODIFS+1
1352 do 5283 j=1,JM0
1353 do 5283 k=1,lmo
1354 RTGOAV(j,k)=RTGOAV(j,k)+RTGO(1,j,k)
1355 5283 continue
1356 endif
1357 C**** RESTRUCTURE THE OCEAN LAYERS AND ELIMINATE SMALL ICE BERGS 203.15
1358 CALL OSTRUC 203.16
1359 if(JDATE.eq.41)then
1360 print *,JDATE,JMONTH,JYEAR
1361 print *,' main after ostruc'
1362 print *,'Z1O'
1363 print 5001,(Z1O(1,j),j=1,JM0)
1364 print *,' T1 ocean'
1365 print 5001,(ODATA(1,j,1),j=1,JM0)
1366 print *,' T2 top ocean'
1367 print 5001,((2.*ODATA(1,j,4)-ODATA(1,j,5)),j=1,JM0)
1368 print *,' T2 ocean'
1369 print 5001,(ODATA(1,j,4),j=1,JM0)
1370 print *,' T3 ocean'
1371 print 5001,(ODATA(1,j,5),j=1,JM0)
1372 endif
1373 ENDIF ! KOCEAN
1374 #endif
1375 #if ( defined ML_2D )
1376 if(TRANSR)then
1377 NODIFS=NODIFS+1
1378 do 5283 j=1,JM0
1379 do 5283 k=1,lmo
1380 RTGOAV(j,k)=RTGOAV(j,k)+RTGO(1,j,k)
1381 5283 continue
1382 endif
1383 #endif
1384
1385 #if ( defined CPL_OCEANCO2 )
1386
1387 C For OCM or 3D ocean model with carbon cycle
1388 #ifdef PREDICTED_GASES
1389 c -------
1390 c 102596
1391 c air co2 mixing ratio goes to ocean carbon model:
1392 c
1393 do j=1,jm
1394 pC_atm(j)=zco2(1,j,1)
1395 & *28.97296245/44.0*1.e-9
1396 !ppb(m) to kg per volume base
1397
1398 atm_co2(j)=pC_atm(j)*1.e6
1399
1400 enddo ! j
1401 c
1402 c -------
1403 #else
1404
1405 do j=1,jm
1406 pC_atm(j)=atm_co2(j)*1.e-6
1407 enddo ! j
1408
1409 #endif
1410
1411 #if ( defined ML_2D )
1412 CB Gary CO2 uptake by the ocean
1413
1414 do j=1,jm
1415 if(NWMGEA(J).gt.0)then
1416 WSAV(J)=WSAV(J)/NWMGEA(J)
1417 NWMGEA(J)=0.
1418 else
1419 WSAV(J)=0.
1420 end if
1421 end do ! j
1422
1423 do j=1,jm
1424 tggary(j)=ODATA(1,j,1)+273.16
1425 wsgary(j)=WSAV(J)
1426 WSAV(j)=0.
1427 arsrf(j)=areaml(j)*(1.-ODATA(1,j,2))
1428 DEPTHML(j)=ZOAV(J)
1429 co24ocnan(j)=co24ocnan(j)+pC_atm(j)*1.e6
1430 enddo ! j
1431 c print *,'CO2 for 2D ocean'
1432 c print *,JYEAR,JMONTH
1433 c print *,'co2=',pC_atm(27)*1.e6,' ws=',wsgary(27)
1434 c print *,'tem=',tggary(27)
1435 c print '12f7.1,/,2(11f7.1,/,),12f7.1',(pC_atm(j)*1.e6,j=1,jm)
1436 c print '12f7.1,/,2(11f7.1,/,),12f7.1',(rco2(j,1),j=1,jm)
1437 c ncallgary=ncallgary+1
1438 ! 10/28/06
1439 ! call carb_mxdlyr_chem(focean)
1440 ! call carb_airsea_flx
1441 !
1442 ! 3D ocean chemistry
1443 call carb_chem_ocmip(focean)
1444 call carb_airsea_flx(dtco2)
1445 ! 3D ocean chemistry
1446 ! 10/28/06
1447 c print *,'FCO2 ncallgary=',ncallgary
1448 c print '12f7.1,/,2(11f7.1,/,),12f7.1',
1449 c & (fluxco2(j)*12.e-15*365.,j=1,jm)
1450 #endif
1451
1452 C For ocean carbon model
1453 c Annual oceanic CO2 uptake
1454 do j=1,jm
1455 OCUPT=OCUPT+fluxco2(j)
1456 enddo
1457 print *,' OCUPT=',OCUPT*12.e-15
1458
1459 #if ( defined CPL_CHEM )
1460 !
1461 ! === Calculate ocean uptake of CO2
1462 ! === once per day:
1463 !
1464 i=1
1465 call chemairmass(p) !update airmass
1466
1467 call chemoceanco2(fluxco2)
1468 !
1469 #endif
1470
1471 do j=1,jm
1472 fluxco2(j)=0.0
1473 enddo
1474 #endif
1475
1476 #if ( defined CPL_TEM )
1477 !#if ( defined CLM )
1478 c print *,'JDATE for TEM=',JDATECLM
1479 do j=1,jm
1480 if(npred4tem(j).gt.0)then
1481 c pred4tem(j)=pred4tem(j)/npred4tem(j)
1482 ewvd4tem(j)=ewvd4tem(j)/npred4tem(j)
1483 pre4tem(J)=pre4tem(J)+pred4tem(j)
1484 endif
1485 c
1486 if(nt2md4tem(j).gt.0)then
1487 t2md4tem(j)=t2md4tem(j)/nt2md4tem(j)
1488 endif
1489 temp4tem(j)=temp4tem(j)+t2md4tem(j)
1490 dtem4tem(JDATECLM,j)=t2md4tem(j)
1491 c
1492 if(nradd4tem(j).gt.0)then
1493 cldd4tem(j)=cldd4tem(j)/ncldd4tem(j)
1494 swtd4tem(j)=swtd4tem(j)/nradd4tem(j)
1495 swsd4tem(j)=swsd4tem(j)/nradd4tem(j)
1496 endif
1497 sws4tem(j)=sws4tem(j)+swsd4tem(j)
1498 c
1499 enddo
1500 c
1501 do j=1,jm
1502 pred4tem(j)=0.0
1503 ewvd4tem(j)=0.0
1504 t2md4tem(j)=0.0
1505 cldd4tem(j)=0.0
1506 swtd4tem(j)=0.0
1507 swsd4tem(j)=0.0
1508 npred4tem(j)=0
1509 ncldd4tem(j)=0
1510 nradd4tem(j)=0
1511 nt2md4tem(j)=0.
1512 enddo
1513 #endif
1514
1515 #if ( defined CPL_OCEANCO2 && defined ML_2D )
1516 C For OCM
1517
1518 ! dtco2=3600.*24.
1519 ! cfkvct=1.0
1520 ! if (JYEAR.ge.1991) then begin
1521 ! if (JYEAR.le.2100) then begin
1522 ! cfkvct=(1.0*(2100-JYEAR)+0.25*(JYEAR-1990))/110.
1523 ! esle
1524 ! cfkvct=0.25
1525 ! endif
1526 ! endif
1527 ! do j=1,jm
1528 ! edzcatr(j)=cfkvct*edzcar(j)
1529 ! enddo
1530 call diffusco2(lmo,jm,dtco2,0.5,edzcart,depthml,focean,
1531 & dzg,dzog,rco2)
1532 call hdocean(rco2,focean,dxv,dyv,DXYP,depthml,edohd,dtco2)
1533 call avegary
1534 CB Gary CO2 uptake by the ocean
1535 #endif
1536
1537 #if ( defined CPL_CHEM) && ( defined CPL_TEM )
1538 C take into accout land uptake form TEM for previous month
1539 do j=1,jm
1540 fluxnep(j)=aduptd+1.e-3*temco2(j)/NDAYMN(mndriver)
1541 c Annual TEM CO2 uptake
1542 TEMUPTANN=TEMUPTANN+fluxnep(j)
1543 enddo
1544 if(jdate.eq.1)then
1545 print *,'Monthly TEM uptake'
1546 c print *,mndriver,adupt,temuptann-temup0
1547 temup0=temuptann
1548 endif
1549 C
1550 c
1551 i=1
1552 call chemairmass(p) !update airmass
1553
1554 call chemtemco2(fluxnep)
1555 C
1556 c
1557 !
1558 #endif
1559
1560 c End of month
1561 if(JDATE.eq.1)then
1562
1563 #if ( defined CPL_CHEM )
1564 !
1565 ! === Calculating monthly averaged mixing ratios:
1566 !
1567 call chemmonth2
1568 c print *,'Atmosphere after chemmonth2'
1569
1570 ! === Calculate and print monthly n-s transport
1571 ! === of cfc11:
1572 !
1573 ! call chemtmp2
1574 !
1575 do nhr = 1,8
1576 sfc3hro3(1:nlat,nhr) = sfc3hro3(1:nlat,nhr)
1577 & /float(cnt3hr(nhr))
1578 #if ( defined CPL_TEM )
1579 o34tem(nhr,1:nlat)=sfc3hro3(1:nlat,nhr)
1580 #endif
1581 end do
1582 cwrite(124)sfc3hro3
1583 cnt3hr(1:8) = 0
1584 sfc3hro3(1:nlat,1:8) = 0.0
1585
1586 ! === Writing rawdata every month:
1587 !
1588 ! call chemprt !closed 032697
1589
1590 call chembudget (p)
1591 print *,' Atmosphre after chembudget mymonth=',mymonth
1592 ! === 09/26/94
1593 ! === Reset year and month index:
1594 !
1595 mymonth = mymonth + 1
1596 if(mymonth.gt.12)then
1597 myyear = myyear +1
1598 ! myyear = min(myyear,nchemyr)
1599 myyear = min(myyear,myyearlast)
1600 mymonth = 1
1601 ! endif ! 27/8/2005
1602
1603 ! === 092295
1604 ! === write rawdata for possible renew run
1605 ! === at end of each month:
1606 ! === at end of each year: 27/8/2005
1607 !
1608 rewind 178
1609 print *,'For chem restart ',myyear,mymonth
1610 write(178)myyear,mymonth,airmass,
1611 & cfc11,cfc110,
1612 & cfc11m,
1613 & cfc11sd,
1614 & cfc12,cfc12m,
1615 & cfc12sd,
1616 & xn2o ,xn2om ,
1617 & xn2osd,
1618 & hfc134a,hfc134am,
1619 & pfc ,pfcm,
1620 & sf6 ,sf6m,
1621 & bcarbon,bcm,
1622 & ocarbon,ocm,
1623 & atomo ,
1624 & o1d ,
1625 & o3 ,o3m ,
1626 & co ,com ,
1627 & zco2 ,zco2m,
1628 & atomh ,
1629 & ho ,
1630 & ho2 ,hoxm ,
1631 & h2o2 ,
1632 & xno ,
1633 & xno2 ,xnoxm,
1634 & xno3 ,
1635 & xn2o5 ,xnoym,
1636 & hno3 ,
1637 & ch4 ,ch4m ,
1638 & ch3 ,
1639 & cho ,
1640 & ch2o ,
1641 & ch3o ,
1642 & ch3o2 ,
1643 & ch3o2h,
1644 & so2 ,so2m ,
1645 & hoso2 ,
1646 & so3 ,
1647 & h2so4 ,h2so4m,
1648 & sviod ,sviodm
1649 rewind 178
1650
1651 endif
1652 #endif
1653
1654 #if ( defined CPL_TEM )
1655 !#if ( defined CLM )
1656 do J=1,JM
1657 c
1658 #ifdef PREDICTED_GASES
1659 co24tem(j)=zco2(1,j,1)
1660 & *28.97296245/44.0*1.e-3
1661 !ppm(m) to kg per volume base
1662 #else
1663 co24tem(j)=atm_co2(j)
1664
1665 #endif
1666 c
1667 enddo
1668 #endif
1669
1670 endif ! end of month
1671 550 continue
1672 C END of DAilY
1673
1674
1675 C CALL CHECKT (11) 203.17
1676 CALL CLOCKS (MNOW) 203.18
1677 MSURF=MSURF+(MLAST-MNOW) 203.19
1678 MLAST=MNOW 203.2
1679 C**** 204.
1680 C**** WRITE INFORMATION ONTO A TAPE EVERY USET HOURS 205.
1681 C**** 206.
1682 IF(USET.LE.0.) GO TO 600 207.
1683 IF(.NOT.EVENT(USET)) GO TO 600 208.
1684 C COMPUTATIONS FOR XXXXXX 209.
1685 WRITE (520) TAU,XXXXXX 210.
1686 CALL CLOCKS (MNOW) 211.
1687 MINC=MLAST-MNOW 212.
1688 MELSE=MELSE+MINC 213.
1689 PERCNT=100.*MELSE/(MSTART-MNOW+1.E-5) 214.
1690 c WRITE (6,910) MINC,MELSE,PERCNT,TAU 215.
1691 C**** 216.
1692 C**** CALL DIAGNOSTIC ROUTINES 217.
1693 C**** 218.
1694 600 IF(MOD(NSTEP-NDYN,NDA4).EQ.0) CALL DIAG4A 219.
1695 C 220.
1696 IF(NDPRNT(1).GE.0) GO TO 610 221.
1697 C**** PRINT CURRENT DIAGNOSTICS (INCLUDING THE INITIAL CONDITIONS) 222.
1698 IF(KDIAG(1).LT.9) CALL DIAG1(NOCLM) 223.
1699 IF(KDIAG(2).LT.9) CALL DIAG2 224.
1700 IF(KDIAG(7).LT.9) CALL DIAG7P 225.
1701 IF(KDIAG(3).LT.9) CALL DIAG3 226.
1702 C 227.
1703 C 228.
1704 IF(KDIAG(4).LT.9) CALL DIAG4 229.
1705 IF(TAU.LE.TAUI+DTHR*(NDYN+.5)) CALL DIAGKN 230.
1706 NDPRNT(1)=NDPRNT(1)+1 231.
1707 C 690 changed to 691 02/21/2003
1708 610 IF(.NOT.EVENT(24.)) GO TO 691 232.
1709 C**** PRINT DIAGNOSTIC TIME AVERAGED QUANTITIES ON NDPRNT-TH DAY OF RUN 233.
1710 DO 620 K=1,13 234.
1711 IF (JDAY.EQ.NDPRNT(K)) GO TO 630 235.
1712 620 CONTINUE 236.
1713 GO TO 640 237.
1714 c 630 WRITE (6,920) 238.
1715 630 continue
1716 IF(KDIAG(1).LT.9) CALL DIAG1(NOCLM) 239.
1717 IF(KDIAG(2).LT.9) CALL DIAG2 240.
1718 IF(KDIAG(7).LT.9) CALL DIAG7P 241.
1719 IF(KDIAG(3).LT.9) CALL DIAG3 242.
1720 C 243.
1721 C 244.
1722 C IF(KDIAG(6).LT.9) CALL DIAG6 245.
1723 IF(KDIAG(4).LT.9) CALL DIAG4 246.
1724 IF(KDIAG(8).LT.9) CALL DIAG8 (IPFLAG) 247.
1725 C**** THINGS TO DO BEFORE ZEROING OUT THE ACCUMULATING ARRAYS 248.
1726 C**** (NORMALLY DONE AT THE END OF A MONTH) 249.
1727 640 DO 650 K=1,13 250.
1728 IF(JDAY.EQ.NDZERO(K)) GO TO 660 251.
1729 650 CONTINUE 252.
1730 GO TO 690 253.
1731 C**** PRINT THE KEY DIAGNOSTICS 254.
1732 660 CONTINUE 255.
1733 CALL DIAGKN 255.
1734 C**** PRINT AND ZERO OUT THE TIMING NUMBERS 256.
1735 CALL CLOCKS (MNOW) 257.
1736 MDIAG=MDIAG+(MLAST-MNOW) 258.
1737 MLAST=MNOW 259.
1738 TOTALT=.01*(MSTART-MNOW) 260.
1739 PDYN=MDYN/TOTALT 261.
1740 PCDNS=MCNDS/TOTALT 262.
1741 PRAD=MRAD/TOTALT 263.
1742 PSURF=MSURF/TOTALT 264.
1743 PDIAG=MDIAG/TOTALT 265.
1744 PELSE=MELSE/TOTALT 266.
1745 DTIME=24.*TOTALT/(60.*(TAU-TAU0)) 267.
1746 c WRITE (6,909) DTIME,PDYN,PCDNS,PRAD,PSURF,PDIAG,PELSE 268.
1747 MDYN=0 269.
1748 MCNDS=0 270.
1749 MRAD=0 271.
1750 MSURF=0 272.
1751 MDIAG=0 273.
1752 MELSE=0 274.
1753 MSTART=MNOW 275.
1754 if(TRANSR)then
1755 do 5284 j=1,JM0
1756 do 5284 k=1,lmo
1757 RTGOAV(j,k)=RTGOAV(j,k)/NODIFS
1758 5284 continue
1759 c print *,'ATM RTGOAV monthly'
1760 c print *,(RTGOAV(J,1),j=1,jm)
1761 endif ! TRANSR
1762 SPGAV=srps(jm+3)/nsrps+PTOP
1763 do 5285 j=1,JM+3
1764 GBUDG(j,38,1)=(srps(j)/nsrps+PTOP)*1013./SPGAV
1765 5285 continue
1766 c do 5287 j=1,JM+3
1767 c GBUDG(j,38,1)=GBUDG(j,37,1)*1013./GBUDG(jm+3,37,1)
1768 c5287 continue
1769 ! print *,'FRMDICE'
1770 ! print '6(1PE12.4)',FRMDICE
1771 ENKE=0.0
1772 ENPT=0.0
1773 do ii=1,4
1774 ENKE=ENKE+SPECA(1,19,ii)
1775 ENPT=ENPT+SPECA(1,20,ii)
1776 enddo
1777 c print *,'ENKE=',ENKE,' ENTP=',ENPT,' ENTT=',ENKE+ENPT
1778 print *,'ENKE=',QTABLE(JMP3,19,10),' ENTP=',QTABLE(JMP3,20,10),
1779 & ' ENTT=',QTABLE(JMP3,19,10)+QTABLE(JMP3,20,10)
1780 print *,'WLMMAX=',WLMMAX,' JWLMMAX=',JWLMMAX
1781 c print *,'AJ(*,37)'
1782 c print *,(AJ(J,37),j=1,jm)
1783 c print *,'AJ(*,28)'
1784 c print *,(AJ(J,28),j=1,jm)
1785 c IF(USEP.LE.0.) GO TO 680 276.
1786 C**** WRITE SELECTED DIAGNOSTICS ONTO A DISK DATA SET FOR PLOTTING 277.
1787 c IF(TAU.LE.TAUI+1080.) GO TO 675 278.
1788 c 670 READ (16) TAUX 279.
1789 c IF(TAU.GT.TAUX+1080.) GO TO 670 280.
1790 675 WRITE (546) AEXP4,JDATE,JMONTH,JYEAR,JDATE0,JMNTH0,JYEAR0, 281.
1791 * GBUDG,QMAPS,QTABLE,INQTAB,J1QT,INQMAP,RTGOAV
1792 print *,'From atm write(546) ',JMNTH0,JYEAR0,' ',JMONTH,JYEAR
1793 c print *,0.1*GBUDG(1,26,2),0.1*GBUDG(1,35,2)
1794 nwr=nwr+1
1795 690 continue
1796 if(JDAY.eq.1)then
1797 rewind 547
1798 write(547)AEXP,nwr
1799 rewind 547
1800 print *,'From atm write(547) ',AEXP,nwr
1801 endif
1802 if(wr25.and.JDAY.eq.1)then
1803 c Write a restart file once a year
1804 print *,'Write a restart file once a year.'
1805
1806 #if ( defined CPL_OCEANCO2 && defined ML_2D )
1807 C Data for possible restart for OCM
1808 write(369)jyear-1,vdfocm
1809 write(369)Hg
1810 write(369)Rco2
1811 rewind 369
1812 #endif
1813
1814 print *,' KDISK0=',KDISK0
1815 CALL RFINAL (IRAND)
1816 REWIND KDISK0
1817 if(TRANSR)then
1818 WRITE(KDISK0) AEXP,TAU,JC,C,RC,KEYNR,U,V,T,P,Q,ODATA,GDATA,
1819 * BLDATA,
1820 * RQT,SRHR,TRHR,(AJ(K,1),K=1,KACC),TAU,TSSFC,CKS,CKN,WMGE,TPRIM2
1821 * ,MRCHT,TRSURF,SRSURF,TLANDW,TSURFW,DWAV0
1822 * ,TG3M,RTGO,STG3,DTG3
1823 else
1824 WRITE(KDISK0) AEXP,TAU,JC,C,RC,KEYNR,U,V,T,P,Q,ODATA,GDATA,
1825 * BLDATA,
1826 * RQT,SRHR,TRHR,(AJ(K,1),K=1,KACC),TAU,TSSFC,CKS,CKN,WMGE,TPRIM2
1827 * ,MRCHT,TRSURF,SRSURF,TLANDW,TSURFW,DWAV0
1828 c print *,' TSURFT'
1829 c print 5001,TSURFT
1830 c print *,' TSURFW'
1831 c print 5001,TSURFW
1832 endif
1833 REWIND KDISK0
1834 KDISK=3-KDISK
1835 KDISK0=500+KDISK
1836 end if
1837 C 690 changed to 691 02/21/2003
1838 680 IF(KCOPY.LE.0) GO TO 691 284.
1839 C**** WRITE A COPY OF THE FINAL RESTART DATA SET ONTO DISK 285.
1840 CALL RFINAL (IRAND) 285.5
1841 print *,' after 680'
1842 print *,' TAU=',TAU,' IRAND=',IRAND
1843 IF(KCOPY.GT.99) GO TO 687 286.
1844 683 READ (KCOPY) TAUX 286.5
1845 IF(TAU.GT.TAUX+3240.) GO TO 683 287.
1846 685 WRITE (KCOPY) TAU,JC,C,RC,KEYNR,U,V,T,P,Q,ODATA,GDATA,BLDATA, 287.5
1847 * RQT,SRHR,TRHR,(AJ(K,1),K=1,KACC),TAU,TSSFC,CKS,CKN 288.
1848 * ,WMGE,TPRIM2,MRCHT,TRSURF,SRSURF,TLANDW,TSURFW,DWAV0
1849 * ,TG3M,RTGO,STG3,DTG3
1850 REWIND KCOPY 288.5
1851 C 690 changed to 691 02/21/2003
1852 GO TO 691 289.
1853 687 KCOPY=KCOPY-100 289.5
1854 GO TO 685 289.6
1855 C**** TIME FOR CALLING DIAGNOSTICS 290.
1856 C 690 changed to 691 02/21/2003
1857 691 CALL CLOCKS (MNOW) 291.
1858 MDIAG=MDIAG+(MLAST-MNOW) 292.
1859 MLAST=MNOW 293.
1860 780 IF(TAU.LE.TAUI+DTHR*(NDYN+.5)) GO TO 120 294.
1861 GO TO 100 295.
1862 C**** 296.
1863 C**** END OF MAIN LOOP 297.
1864 C**** 298.
1865 C**** RUN TERMINATED BECAUSE SENSE SWITCH 6 WAS TURNED ON 299.
1866 800 WRITE (6,904) 300.
1867 IF(EVENT(TAUT)) GO TO 820 301.
1868 CALL RFINAL (IRAND) 302.
1869 print *,' after 800'
1870 print *,' TAU=',TAU,' IRAND=',IRAND
1871 if(wr25) then
1872 REWIND KDISK0 303.
1873 WRITE(KDISK0) AEXP,TAU,JC,C,RC,KEYNR,U,V,T,P,Q,ODATA,GDATA,BLDATA, 304.
1874 * RQT,SRHR,TRHR,(AJ(K,1),K=1,KACC),TAU,TSSFC,CKS,CKN,WMGE,TPRIM2 305.
1875 * ,MRCHT,TRSURF,SRSURF,TLANDW,TSURFW,DWAV0
1876 end if
1877 C WRITE (6,908) 306.
1878 C**** RUN TERMINATED BECAUSE IT REACHED TAUE (OR SS6 WAS TURNED ON) 307.
1879 c 820 WRITE (6,905) TAU,IDAY,TOFDAY 308.
1880 820 continue
1881
1882 #if ( defined OCEAN_3D || defined ML_2D )
1883 C DTATM time step of atm model in hours
1884 C precip and evap in mm/day or kg/m**2/day
1885 do j=1,jm0
1886 #if ( defined OCEAN_3D && defined CPL_OCEANCO2 )
1887 ncallatm=ncallatm+1
1888 ! 020107
1889 ! co24ocean(j)=pC_atm(j)*1.e6
1890 co24ocean(j)=atm_co2(j)
1891 ! 020107
1892 co24ocnan(j)=co24ocnan(j)+co24ocean(j)
1893 #endif
1894 #ifdef ML_2D
1895 cjrs block only MD_2D
1896 rseaice(j)=ODATA(1,J,2)
1897 #endif
1898 tauu(j)=tauu(j)/(NSURF*DTATM)
1899 tauv(j)=tauv(j)/(NSURF*DTATM)
1900 tempr(j)=tempr(j)/DTATM
1901 precip(j)=precip(j)/(DTATM/24.)
1902 fland=FDATA(1,J,2)
1903 if (fland.lt.1.0)then
1904 precip(j)=precip(j)*(1.-fland*prlnd2total(j,mndriver))
1905 & /(1.-fland)
1906 endif
1907 Cjrs surfpr(j)=surfpr(j)/(DTATM/24.)
1908 c
1909 if(naveo(j).gt.0)then
1910 hfluxo(j)=NSURF*hfluxo(j)/(NDYN*DT*naveo(j))
1911 dhfodtg(j)=dhfodtg(j)/naveo(j)
1912 dhfodtgeq(j)=dhfodtgeq(j)/naveo(j)
1913 evao(j)=NSURF*evao(j)/(NDYN*DT*naveo(j))
1914 devodtg(j)=devodtg(j)/naveo(j)
1915 devodtgeq(j)=devodtgeq(j)/naveo(j)
1916 C From mm/sec to mm/day
1917 evao(j)=24.*3600.*evao(j)
1918 devodtg(j)=24.*3600.*devodtg(j)
1919 devodtgeq(j)=24.*3600.*devodtgeq(j)
1920 endif
1921 C
1922 if(navei(j).gt.0)then
1923 hfluxi(j)=NSURF*hfluxi(j)/(NDYN*DT*navei(j))
1924 dhfidtg(j)=dhfidtg(j)/navei(j)
1925 dhfidtgeq(j)=dhfidtgeq(j)/navei(j)
1926 evai(j)=NSURF*evai(j)/(NDYN*DT*navei(j))
1927 devidtg(j)=devidtg(j)/navei(j)
1928 devidtgeq(j)=devidtgeq(j)/navei(j)
1929 C From mm/sec to mm/day
1930 evai(j)=24.*3600.*evai(j)
1931 devidtg(j)=24.*3600.*devidtg(j)
1932 devidtgeq(j)=24.*3600.*devidtgeq(j)
1933 endif
1934 C
1935 if(navrad(j).gt.0)then
1936 solarinc_ice(j)=solarinc_ice(j)/navrad(j)
1937 solarnet_ice(j)=solarnet_ice(j)/navrad(j)
1938 endif
1939 if(navrado(j).gt.0)then
1940 solarinc_ocean(j)=solarinc_ocean(j)/navrado(j)
1941 solarnet_ocean(j)=solarnet_ocean(j)/navrado(j)
1942 endif
1943 c Runoff is a flux of water from land in mm/day
1944 c not for m**2
1945 cjrs change runoff to new name arunoff
1946 arunoff(j)=arunoff(j)/(DTATM/24.)*FDATA(1,j,2)
1947 & *DXYP(J)
1948 if(NWMGEA(J).gt.0)then
1949 wsocean(J)=WSAV(J)/NWMGEA(J)
1950 NWMGEA(J)=0.
1951 else
1952 wsocean(J)=0.
1953 end if
1954 WSAV(J)=0.
1955 c dhdtav(j)=dhdtav(j)+dhfdtg(j)
1956 c devdtav(j)=devdtav(j)+devdtg(j)
1957 c
1958 c
1959 ps4ocean(j)=ps4ocean(j)/DTATM
1960 do l=1,lm
1961 qyz4ocean(j,l)=qyz4ocean(j,l)/DTATM
1962 tyz4ocean(j,l)=tyz4ocean(j,l)/DTATM
1963 enddo
1964 c
1965 c
1966 end do ! j
1967 rungl=0.0
1968 runn=0.0
1969 runt=0.0
1970 runs=0.0
1971 SLAND=0.0
1972 CLAT=20.*TWOPI/360.
1973 do j=1,jm
1974 SLAND=SLAND+FDATA(1,j,2)*DXYP(J)
1975 c jrs runoff->arunoff
1976 rungl=rungl+arunoff(j)
1977 if(LAT(J).lt.-CLAT)then
1978 runs=runs+arunoff(j)
1979 else if(LAT(J).lt.CLAT)then
1980 runt=runt+arunoff(j)
1981 else
1982 runn=runn+arunoff(j)
1983 endif
1984 enddo
1985 c print *,'RUNOFF TOFDAY=',TOFDAY
1986 c print *,rungl/SLAND,rungl,runs,runt,runn
1987 c nmonth=JMNTH0
1988 #ifdef ML_2D
1989 c jrs only ML_2D
1990 nmonth=AMONTH(mndriver)
1991 #endif
1992 jdatefl=jdate-1
1993 c if(JMONTH.ne.JMNTH0)jdatefl=NDAYMN(mndriver)
1994 c if(JMONTH.ne.JMNTH0)then
1995 c print *,'OCEAN FLUXS'
1996 c print *,JMONTH,JMNTH0,JDATE
1997 c print *,'Month=',AMONTH(mndriver),' day=',jdatefl
1998 c print *,' TAIR_OCEAN'
1999 c print *,(tairo(j),j=1,jm0)
2000 c print *,' TAIR_ICE'
2001 c print *,(tairi(j),j=1,jm0)
2002 c print *,' TAUU'
2003 c print *,(tauu(j),j=1,jm0)
2004 c print *,' TAUV'
2005 c print *,(tauv(j),j=1,jm0)
2006 c print *,' EVA-E OCEAN'
2007 c print *,(evao(j),j=1,jm0)
2008 c print *,' EVA-I OCEAN'
2009 c print *,(evai(j),j=1,jm0)
2010 c print *,' P-E ICE'
2011 c print *,(pmei(j),j=1,jm0)
2012 c print *,' HEAT FLUX OCEAN'
2013 c print *,(hfluxo(j),j=1,jm0)
2014 c print *,' DHEATFLUX_OCEAN/DTG'
2015 c print *,(dhfodtg(j),j=1,jm0)
2016 c print *,' EVA_OCEAN/DTG'
2017 c print *,(devodtg(j),j=1,jm0)
2018 c print *,' EVA_ICE/DTG'
2019 c print *,(devidtg(j),j=1,jm0)
2020 c print *,' HEAT FLUX ICE'
2021 c print *,(hfluxi(j),j=1,jm0)
2022 c print *,' DHEATFLUX_ICE/DTG'
2023 c print *,(dhfidtg(j),j=1,jm0)
2024 c print *,' DLH_OCEAN/DTG'
2025 c print *,(devidtg(j),j=1,jm0)
2026 c print *,'PS4OCEAN'
2027 c print *,(ps4ocean(j),j=1,jm0)
2028 c print *,'QYZ4OCEAN'
2029 c do l=1,lm
2030 c print *,(qyz4ocean(j,l),j=1,jm0)
2031 c enddo
2032 c endif
2033 c go to 587
2034 c write(893),nmonth,jdatefl,tempr,tauu,tauv,precip,evao,
2035 c & evai,hfluxo,dhfodtg,devodtg,hfluxi,dhfidtg,devidtg,
2036 c & solarinc_ice,solarnet_ice,rseaice
2037 #ifdef ML_2D
2038 c jrs only ML_2D
2039 do j=1,jm
2040 osst(j)=ODATA(1,j,1)
2041 aoice(j)=ODATA(1,j,3)
2042 foice(j)=ODATA(1,j,2)
2043 snowice(j)=GDATA(1,j,1)
2044 tice1(j)=GDATA(1,j,3)
2045 tice2(j)=GDATA(1,j,7)
2046 enddo
2047 #endif
2048 c write (894),nmonth,jdatefl,osst,aoice,foice,snowice,tice1,tice2
2049 587 continue
2050 #endif
2051
2052 C
2053 c TAU for coupler
2054 c
2055 TAUATM=TAU
2056 MONTHATM=JMONTH
2057 JDATEATM=JDATE
2058 JYEARATM=JYEAR
2059 C
2060 #ifdef ML_2D
2061 Cjrs change this block to only ML_2D
2062 IDAYM=IDAY
2063 JDAYM=JDAY
2064 JDATEM=JDATE
2065 JMONTHM=JMONTH
2066 JYEARM=JYEAR
2067 TAUML=TAU
2068 TOFDAYML=TOFDAY
2069 #endif
2070 C
2071 if(JDAY.ne.JDAYLAST)then
2072 c print *,'co24ocean=',co24ocean(jm/2)
2073 if(JDATEPR.ne.0)then
2074 c WRITE (6,905) TOFDAY,JDATE,JMONTH,JYEAR
2075 WRITE (6,905) TOFDAYPR,JDATEPR,JMONTHPR,JYEARPR
2076 endif
2077 cjrs print *,'ncallclm=',ncallclm
2078 JDAYLAST=JDAY
2079 c if(ncallclm.gt.6) stop
2080 c stop
2081 endif
2082 return
2083 C CALL ENQJOB 309.
2084 C CALL ENQJOB 310.
2085 IF(KSS6.EQ.1) STOP 12 310.1
2086 IF(IPFLAG.EQ.0) STOP 13 311.
2087 STOP 1 312.
2088 C**** 313.
2089 901 FORMAT ('0CLIMATE MODEL STARTED UP',14X,'DAY',I5,', HR',F6.2,I6, 314.
2090 * A5,I27,I7,F7.1,' TAU',F9.2) 315.
2091 902 FORMAT (' DYNAMIC TERMS INTEGRATED, MRCH=',I1,4X,'DAY',I5, 316.
2092 * ', HR',F6.2,I6,A5,2I7,F7.1,23X,'TAU',F9.2) 317.
2093 903 FORMAT (' SOURCE TERMS INTEGRATED',64X,2I7,F7.1) 318.
2094 904 FORMAT ('0SENSE SWITCH 6 HAS BEEN TURNED ON.') 319.
2095 905 FORMAT (/1(1X,33('****')/)/ 320.
2096 c * ' PROGRAM HAS TERMINATED NORMALLY. TAU,IDAY,TOFDAY=',F9.2, 321.
2097 * ' ATM HAS TERMINATED NORMALLY. AT ',F9.2,I3,A5,i5,
2098 * /1(1X,33('****')/)) 322.
2099 906 FORMAT (' OUTPUT RECORD WRITTEN ON UNIT',I3,55X,2I7,F7.1, 323.
2100 * ' TAU',F9.2,' ON ',A4) 324.
2101 908 FORMAT (' OUTPUT RECORD WRITTEN ON UNIT',I3,79X,'TAU',F9.2, 325.
2102 * ' ON ',A4) 326.
2103 909 FORMAT (/'0TIME',F7.2,'(MINUTES) DYNAMICS',F5.1, 327.
2104 * ' CONDENSATION',F5.1,' RADIATION',F5.1,' SURFACE', 328.
2105 * F5.1,' DIAGNOSTICS',F5.1,' OTHER',F5.1//) 329.
2106 910 FORMAT (' INFORMATION WRITTEN ON UNIT 20',57X,2I7,F7.1, 330.
2107 * ' TAU',F9.2,' ON TAPE') 331.
2108 920 FORMAT ('1'/64(1X/)) 332.
2109 END 333.

  ViewVC Help
Powered by ViewVC 1.1.22