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

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

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


Revision 1.3 - (hide annotations) (download)
Wed Sep 13 15:23:04 2006 UTC (18 years, 10 months ago) by jscott
Branch: MAIN
Changes since 1.2: +1 -1 lines
comment debugging print statement

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

  ViewVC Help
Powered by ViewVC 1.1.22