/[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.8 - (hide annotations) (download)
Fri Aug 24 19:26:15 2007 UTC (17 years, 11 months ago) by jscott
Branch: MAIN
Changes since 1.7: +6 -3 lines
have atm give CO2 to ocn even if ocn carbon off

1 jscott 1.5 c source 2007 sokolov users 76203 Apr 25 15:29 atmosphere.F
2 jscott 1.1
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 jscott 1.4 # include "CLM.h"
113     # if ( defined CPL_TEM )
114 jscott 1.1 C For TEM
115 jscott 1.4 # include "TEM.h"
116     # endif
117 jscott 1.1 #endif
118 jscott 1.4
119 jscott 1.1 #if ( defined CPL_OCEANCO2 )
120 jscott 1.6 #include "OCM.h"
121 jscott 1.1 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 jscott 1.4 common/Garydiff/depthml(jm0),edzcar(jm0),dzg(lmo),dzog(lmo-1),
125 jscott 1.1 &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 jscott 1.4 common /Garyvlog/odifcarbon,ocarcont
130     common /Garykvct/cfkvct,edzcart(jm0)
131 jscott 1.1 # endif
132     #endif
133    
134 jscott 1.4 INTEGER dtatm, mndriver !routine arguments jrs
135 jscott 1.1
136     #if ( defined OCEAN_3D || defined ML_2D )
137 jscott 1.2 #include "AGRID.h"
138 jscott 1.4 Cjrs elimated COM file/moved elsewhere#include "HRD4OCN.COM"
139 jscott 1.1 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 jscott 1.4 ! common/TSUR/TSURFC(JM0,0:13),TSURFT(JM0),TSURFD(JM0),DTSURF(JM0)
146     #include "TSRF.COM"
147     real TSURFW(JM0),TLANDW(JM0)
148 jscott 1.1 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 jscott 1.4 dtco2=3600.*24.
259 jscott 1.1 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 jscott 1.4 ! 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 jscott 1.1 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 jscott 1.7 cjrs changed to dfloat 8/2/07
369     TAU=DFLOAT(ITAU)/XINT 31.
370 jscott 1.1 IDAY=1+ITAU/I24 32.
371     TOFDAY=(ITAU-(IDAY-1)*I24)/XINT 33.
372 jscott 1.4 ! if(ISTART.eq.2.or.ISTRT1.eq.0.and..not.CONTRR)then
373     ! do 458 j=1,JM
374     ! TSURFD(j)=0.
375     ! TSURFT(j)=0.
376     ! 458 continue
377     ! endif
378 jscott 1.1 if(JDATE.eq.100)then
379     print *,JDATE,JMONTH,JYEAR
380     print *,' main before daily0'
381     print *,' T1 ocean'
382     print 5001,(ODATA(1,j,1),j=1,JM0)
383     print *,' T2 ocean'
384     print 5001,(ODATA(1,j,4),j=1,JM0)
385     print *,' T3 ocean'
386     print 5001,(ODATA(1,j,5),j=1,JM0)
387     print *,' sea ice'
388     print 5002,(ODATA(1,j,2),j=1,JM0)
389     endif
390     CALL DAILY_NEW0 34.
391     print *,' Main after DAILYNEW0 JYEAR=',JYEAR
392 jscott 1.4 print *,"DT2MGL"
393     print *,DT2MGL
394     print *,"DT2MLD"
395     print *,DT2MLD
396 jscott 1.1 #if( !defined OCEAN_3D&& !defined ML_2D )
397     CALL DAILY_OCEAN
398     print *,' AFTER DAILY_OCEAN IDAY=',IDAY,' IYEAR=',IYEAR
399     print *,' JYEAR=',JYEAR,' JDAY=',JDAY
400     print *,' JDATE=',JDATE,' JMONTH=',JMONTH
401     #endif
402     if(JDATE.eq.100)then
403     print *,JDATE,JMONTH,JYEAR
404     print *,' main after daily0'
405     print *,' T1 ocean'
406     print 5001,(ODATA(1,j,1),j=1,JM0)
407     print *,' T2 ocean'
408     print 5001,(ODATA(1,j,4),j=1,JM0)
409     print *,' T3 ocean'
410     print 5001,(ODATA(1,j,5),j=1,JM0)
411     print *,' sea ice'
412     print 5002,(ODATA(1,j,2),j=1,JM0)
413     endif
414     c print *,' Z12O'
415     c print *,(Z12O(1,j),j=1,JM0)
416     c99 CONTINUE 34.993
417     CALL CLOCKS (MLAST) 35.
418     MINC=MNOW-MLAST 36.
419     MELSE=MELSE+MINC 37.
420     PERCNT=100.*MELSE/(MSTART-MLAST+1.E-5) 38.
421     c WRITE (6,901) IDAY,TOFDAY,JDATE,JMONTH,MINC,MELSE,PERCNT,TAU 39.
422     DOPK=1. 40.
423     MODD5K=1000 41.
424     IF(TAU.GE.TAUE) GO TO 820 42.
425     HPRNT=.TRUE.
426     HPRNT=.FALSE.
427     JPR=7
428     JPR=1
429     JPR=14
430     LPR=1
431     print *,' MAIN MRCHT=',MRCHT
432     c
433     c print *,' TAUE=',TAUE
434     TAUEM=TAUE
435     FIRST=.FALSE.
436     C
437     c TAU for coupler
438     c
439     TAUATM=TAU
440     MONTHATM=JMONTH
441     JDATEATM=JDATE
442     JYEARATM=JYEAR
443     C
444 jscott 1.4 #if ( defined CPL_OCEANCO2 )
445     do j=1,jm
446     fluxco2(j)=0.0
447     enddo
448     #endif
449 jscott 1.1 #if ( defined CPL_CHEM) && ( defined CPL_TEM )
450     C For TEM
451     if(ISTRT1.eq.0) then
452     c New run
453     c Reading from flin_nep
454     read(537)adupt,temco2
455     else
456     c Restart of the run
457     c Reading from last_nep
458     read(367)adupt,temco2
459     ! & ,temch4,temn2o
460     rewind 367
461     endif
462     !
463     ! adupt= 1.459814341652516
464     ! adupt= 0.9078891180588442
465     ! adupt= 0.25
466     ! adupt= -0.1123070421398009
467     !
468 jscott 1.4 ! adupt= adupt+0.9 ! for vs23
469    
470 jscott 1.1 aduptd=adupt/(365.*JM)
471     temnepgl=0.0
472     do j=1,jm
473     temnepgl=temnepgl+temco2(j)
474     enddo
475     print *,'ADNEP=',adupt
476     print *,'Initial NEP=',adupt+temnepgl*1.e-3
477     temup0=0.0
478     #endif
479     #if (!defined OCEAN_3D && !defined ML_2D)
480     if(TRANSR)then
481     if(LMO.eq.11) then
482     call ODIFS
483     elseif(LMO.eq.12) then
484     call ODIFS12
485     else
486 jscott 1.4 Print *,' Wrong LMO',LMO
487 jscott 1.1 stop
488     endif
489     endif
490     #endif
491 jscott 1.4 !#if (defined PREDICTED_GASES)
492 jscott 1.1 #if (defined CPL_TEM || defined CPL_OCEANCO2 )
493     if(OBSFOR) then
494     call obsco2(iyear,imontha)
495     mnobco2=imonth
496     endif
497     #endif
498 jscott 1.4 !#endif
499     CJRS removed below from ocean_3d
500 jscott 1.1 #ifdef ML_2D
501     do j=1,jm
502     do i=1,io
503     CLAND4OCEAN(i,j)= C3LAND(I,J)
504     enddo
505     enddo
506     IDAYM=IDAY
507     JDAYM=JDAY
508     JDATEM=JDATE
509     JMONTHM=JMONTH
510     JYEARM=JYEAR
511     TAUML=TAU
512     TOFDAYML=TOFDAY
513     KOCEANM=KOCEAN
514     #endif
515     #if (defined OCEAN_3D || defined ML_2D)
516     do l=1,lm
517     sigfl(l)=sig(l)
518     enddo
519     print *,'SIGFL'
520     print *,sigfl
521     #endif
522     #if ( defined CPL_CHEM )
523     do j=1,jm
524     atm_co2(j)=zco2(1,j,1)
525     & *28.97296245/44.0*1.e-9
526     !ppb(m) to kg per volume base
527     & *1.e6
528     enddo
529     #else
530     if(.not.OBSFOR) then
531     CFF=1.
532     if (CO2.gt.0.0)CFF=CO2
533     do j=1,jm
534     atm_co2(j)=CFF*GHGBGR(1)
535     enddo
536     endif
537     #endif
538 jscott 1.4 #if (defined CPL_TEM || defined CPL_OCEANCO2 )
539 jscott 1.1 print *,'ATM_CO2'
540     print *,atm_co2
541 jscott 1.4 #endif
542 jscott 1.1 JDAYLAST=-1
543     ncallclm=0
544     NOCLM=.true.
545     #if ( defined CLM )
546     NOCLM=.false.
547     #endif
548     print *,' atmosphere DTATM=',DTATM
549 jscott 1.4 print *,' It is running'
550 jscott 1.1 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 jscott 1.4 cjrs change var name to arunoff
585 jscott 1.1 arunoff(j)=0.
586     solarinc_ice(j)=0.
587     solarnet_ice(j)=0.
588     solarinc_ocean(j)=0.
589     solarnet_ocean(j)=0.
590 jscott 1.5 Cjrs not used anymore (?) surfpr(j)=0.
591 jscott 1.1 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 jscott 1.5 #ifdef OCEAN_3D
605     C get data from atm-ocean common block
606     do j=1,jm0
607     ODATA(1,j,1)=mmsst(j)
608     ODATA(1,j,2)=mmfice(j)
609     GDATA(1,j,3)=mmtice(j)
610     GDATA(1,j,1)=mmsnowm(j)
611     ODATA(1,j,3)=mmicem(j)
612     GDATA(1,j,7)=0.5*(mmtice2(j)+mmtice1(j))
613     # ifdef CPL_OCEANCO2
614     fluxco2(j)=fluxco2(j) + dtatm*3600.*mmco2flux(j)
615     # endif
616 jscott 1.4 enddo
617 jscott 1.5 #endif
618 jscott 1.1 WLMMAX=0.0
619     C
620     100 IF(.NOT.EVENT(TAUT)) GO TO 200 46.
621     c HPRNT=TAU.ge.17520.00
622     NSTEP1=NSTEP 46.5
623     C**** WRITE RESTART INFORMATION ONTO DISK 47.
624     120 CALL RFINAL (IRAND) 48.
625     IF(NSTEP.EQ.NSTEP2) GO TO 116 48.3
626     DO 115 K=1,22 48.5
627     DO 115 J=1,JM 48.6
628     DO 115 I=2,IO 48.7
629     115 GDATA(I,J,K)=GDATA(1,J,K) 48.9
630     116 CONTINUE 48.91
631     if(wr25.and.ISTART.eq.2)then
632     print *,' main write'
633     print *,' T1 ocean'
634     print 5001,(ODATA(1,j,1),j=1,JM0)
635     print *,' T2 ocean'
636     print 5001,(ODATA(1,j,4),j=1,JM0)
637     print *,' T3 ocean'
638     print 5001,(ODATA(1,j,5),j=1,JM0)
639     REWIND KDISK0 49.
640     if(TRANSR)then
641     WRITE(KDISK0) AEXP,TAU,JC,C,RC,KEYNR,U,V,T,P,Q,ODATA,GDATA,BLDATA, 50.
642     * RQT,SRHR,TRHR,(AJ(K,1),K=1,KACC),TAU,TSSFC,CKS,CKN,WMGE,TPRIM2 51.
643 jscott 1.4 * ,MRCHT,TRSURF,SRSURF,TLANDW,TSURFW,DWAV0
644 jscott 1.1 * ,TG3M,RTGO,STG3,DTG3
645     print *,' STG3'
646     print 5001,(STG3(1,j),j=1,JM0)
647     print *,' DTG3/356'
648     print 5001,(DTG3(1,j)/365.,j=1,JM0)
649     print *,' RTGO'
650     print 5001,((RTGO(1,j,k),j=1,JM0),k=1,lmo)
651     5001 format(23f6.1)
652     5002 format(23f6.3)
653     else
654     WRITE(KDISK0) AEXP,TAU,JC,C,RC,KEYNR,U,V,T,P,Q,ODATA,GDATA,BLDATA, 50.
655     * RQT,SRHR,TRHR,(AJ(K,1),K=1,KACC),TAU,TSSFC,CKS,CKN,WMGE,TPRIM2 51.
656 jscott 1.4 * ,MRCHT,TRSURF,SRSURF,TLANDW,TSURFW,DWAV0
657 jscott 1.1 endif
658     REWIND KDISK0 52.
659     end if ! ISTART.eq.2
660     CALL CLOCKS (MNOW) 53.
661     MINC=MLAST-MNOW 54.
662     MELSE=MELSE+MINC 55.
663     PERCNT=100.*MELSE/(MSTART-MNOW+1.E-5) 56.
664     MLAST=MNOW 59.
665     C**** TEST FOR TERMINATION OF RUN 60.
666     200 READ (503,END=210) LABSSW 61.
667     c HPRNT=TAU.gt.45.0.and.TAU.lt.60.0
668     c HPRNT=TAU.gt.470.0.and.TAU.lt.550.0
669     NCOMP=0
670     210 REWIND 503 61.1
671     IF(LABSSW.EQ.LABEL1) KSS6=1 61.2
672 jscott 1.4 IF(KSS6.EQ.1) GO TO 800 62.
673 jscott 1.1 IF(TAU+.06125.GE.TAUE) GO TO 820 63.
674     JDAY00=JDAY
675     C**** IF TIME TO ZERO OUT DIAGNOSTIC ACCUMULATING ARRAYS, DO SO 64.
676     C**** (NORMALLY DONE AT THE BEGINNING OF A MONTH) 65.
677     IF (TAU.EQ.TAUI) GO TO 260 66.
678     IF(.NOT.EVENT(24.)) GO TO 290 67.
679     DO 250 K=1,13 68.
680     IF(JDAY.EQ.NDZERO(K)) GO TO 260 69.
681     250 CONTINUE 70.
682     GO TO 290 71.
683     260 CONTINUE
684     TAU0=TAU 72.
685     IDAY0=IDAY 73.
686     TOFDY0=TOFDAY 74.
687     JDATE0=JDATE 75.
688     JMNTH0=JMONTH 76.
689     JYEAR0=JYEAR 77.
690     DO 270 I=1,12 78.
691     270 IDACC(I)=0 79.
692     NODIFS=0
693     nsrps=0.
694     DO 280 K=1,KACC 80.
695     280 AJ(K,1)=0. 81.
696     do 5280 j=1,JM
697     tprmg(j)=0.
698     ntprmg(j)=0.
699     Z1OAV(j)=0.
700     NZ1OAV(j)=0
701     NFLIO(J)=0
702     FLIO(J)=0.
703     NCLR(J)=0
704     do 5282 k=1,lmo
705     RTGOAV(j,k)=0.
706     5282 continue
707     do 5281 n=1,12
708     AJCLR(J,n)=0.
709     BJCLR(J,n)=0.
710     CJCLR(J,n)=0.
711     5281 continue
712     5280 continue
713     do 5286 j=1,jm+3
714     srps(j)=0.
715     5286 continue
716    
717     c#if ( defined CPL_OCEANCO2 && defined ML_2D )
718     c call zerogary(ncallgary)
719     c call zerogary
720     c#endif
721    
722     CALL DIAG9A (1) 82.
723    
724     #if ( defined CPL_NEM )
725     C For Emission
726     c === 031097
727     DO 5655 MONTH=1,12
728     IF(JDAY.LE.JDOFM(MONTH+1)) GO TO 5656
729     5655 CONTINUE
730     5656 MNHTEM=MONTH-1
731     if(MNHTEM.eq.0)MNHTEM=12
732     do J=1,JM
733     c ECH4CHIEN(J)=temch4(J)/NDAYMN(MNHTEM)/1000.
734     c EN2OCHIEN(J)=temn2o(J)/NDAYMN(MNHTEM)/1000.
735    
736     ECH4CHIEN(J)=temch4(J)/NDAYMN(MNHTEM)
737     EN2OCHIEN(J)=temn2o(J)/NDAYMN(MNHTEM)
738    
739     enddo
740    
741     C For Emission
742     #endif
743    
744     290 CONTINUE 83.
745     C**** 84.
746     C**** INTEGRATE DYNAMIC TERMS 85.
747     if(HPRNT)then
748     print *,' main before comp1 1',' TAU=',TAU,JDATE,JMONTH
749     #include "PRNT.COM"
750     print *,' GDATA(1,7,5)=',GDATA(1,7,5),' GDATA(1,7,6)='
751     * ,GDATA(1,7,6)
752     endif
753     C**** 86.
754     MODD5D=MOD(NSTEP,NDA5D) 87.
755     IF(MODD5D.EQ.0) CALL DIAG5A (2,0) 88.
756     c IF(NDYNO.EQ.1) then
757     c print *,NDYNO
758     c print *,jm,im,LMT3P1
759     c endif
760     DO 310 J=1,JM 89.
761     DO 300 L=1,LMT3P1 90.
762     DO 300 I=1,IM 91.
763     UX(I,J,L)=U(I,J,L) 92.
764     300 UT(I,J,L)=U(I,J,L) 93.
765     DO 310 L=1,LM 94.
766     DO 310 I=1,IM 95.
767     310 QT(I,J,L)=Q(I,J,L) 96.
768     if(HPRNT)then
769     print *,' main before comp1 2',' TAU=',TAU
770     #include "PRNT.COM"
771     print *,' GDATA(1,7,5)=',GDATA(1,7,5),' GDATA(1,7,6)='
772     * ,GDATA(1,7,6)
773     endif
774     C**** INITIAL FORWARD STEP, QX = Q + .667*DT*F(Q) 97.
775     NS=0 98.
776     MRCH=0 99.
777     CALL COMP1 (UX,VX,TX,PX,Q,U,V,T,P,Q,DTFS,NS) 100.
778     if(HPRNT)then
779     print *,' main after comp1',' TAU=',TAU,' MRCH=',MRCH
780     #include "PRNT.COM"
781     endif
782     C
783     IF(NDYNO.EQ.1) GO TO 320 101.
784     C**** INITIAL BACKWARD STEP IS ODD, QT = QT + DT*F(QX) 102.
785     MRCH=-1 103.
786     CALL COMP1 (UT,VT,TT,PT,QT,UX,VX,TX,PX,Q,DT,NS) 104.
787     if(HPRNT)then
788     print *,' main after comp1',' TAU=',TAU,' MRCH=',MRCH
789     #include "PRNT.COM"
790     endif
791     C
792     GO TO 360 105.
793     C**** INITIAL BACKWARD STEP IS EVEN, Q = Q + DT*F(QX) 106.
794     320 NS=1 107.
795     MODD5K=MOD(NSTEP+NS-NDYN+NDA5K,NDA5K) 108.
796     MRCH=1 109.
797     CALL COMP1 (U,V,T,P,Q,UX,VX,TX,PX,QT,DT,NS) 110.
798     if(HPRNT)then
799     print *,' main after comp1',' TAU=',TAU,' MRCH=',MRCH
800     #include "PRNT.COM"
801     endif
802     C
803     CD DIAGA SHOULD BE CALLED HERE BUT THEN ARRAYS MUST BE CHANGED 111.
804     c
805     c
806     C**** ODD LEAP FROG STEP, QT = QT + 2*DT*F(Q) 112.
807     340 MRCH=-2 113.
808     CALL COMP1 (UT,VT,TT,PT,QT,U,V,T,P,Q,DTLF,NS) 114.
809     if(HPRNT)then
810     print *,' main after comp1',' TAU=',TAU,' MRCH=',MRCH
811     #include "PRNT.COM"
812     endif
813     C
814     C**** EVEN LEAP FROG STEP, Q = Q + 2*DT*F(QT) 115.
815     360 NS=NS+2 116.
816     MODD5K=MOD(NSTEP+NS-NDYN+NDA5K,NDA5K) 117.
817     MRCH=2 118.
818     CALL COMP1 (U,V,T,P,Q,UT,VT,TT,PT,QT,DTLF,NS) 119.
819     if(HPRNT)then
820     print *,' main after comp1',' TAU=',TAU,' MRCH=',MRCH
821     #include "PRNT.COM"
822     endif
823     C
824     IF(NS.LT.NDYN) GO TO 340 122.
825     c IF(MOD(NSTEP+NS-NDYN+NDAA,NDAA).LT.MRCH) THEN 120.
826     CALL DIAGA (UT,VT,TT,PT,QT,NOCLM) 121.
827     if(HPRNT)then
828     print *,' main after DIAGA',' TAU=',TAU,' MRCH=',MRCH
829     #include "PRNT.COM"
830     endif
831     c ENDIF
832     c IF(NS.LT.NDYN) GO TO 340 122.
833     DOPK=1. 123.
834     CALL CLOCKS (MNOW) 124.
835     MINC=MLAST-MNOW 125.
836     MDYN=MDYN+MINC 126.
837     MLAST=MNOW 127.
838     PERCNT=100.*MDYN/(MSTART-MNOW+1.E-5) 128.
839     C 130.
840     CALL DIAG9A (2) 131.
841     C**** 133.
842     C**** INTEGRATE SOURCE TERMS 134.
843     C**** 135.
844     #if (!defined PREDICTED_GASES)
845     #if (defined CPL_TEM || defined CPL_OCEANCO2 )
846     if(OBSFOR) then
847     if(JMONTH.ne.AMONTH(mnobco2)) then
848     mnobco2=mnobco2+1
849     if(mnobco2.eq.13)mnobco2=1
850     call obsco2(jyear,mnobco2)
851     endif
852     endif
853     #endif
854     #endif
855     C
856     C
857     MODRD=MOD(NSTEP,NRAD) 136.
858     MODD5S=MOD(NSTEP,NDA5S) 137.
859     IF(MODD5S.EQ.0) IDACC(8)=IDACC(8)+1 138.
860     C 139.
861     C**** CONDENSTATION, SUPER SATURATION AND MOIST CONVECTION 140.
862     if(HPRNT)then
863     print *,' main before conds',' TAU=',TAU
864     #include "PRNT.COM"
865     endif
866    
867    
868     #if ( defined CPL_CHEM )
869     !
870     ! ===== Calculate airmass at grid point
871     ! ===== Chien Wang, 092395
872     !
873    
874     call chemairmass(p4chem0)
875    
876     !
877     ! === Calculating total mass of tracers
878     ! === with long residence times:
879     !
880     call chemmass1(cfc11, cfc11mass)
881     call chemmass1(cfc12, cfc12mass)
882     call chemmass1(xn2o, xn2omass)
883     call chemmass1(zco2, zco2mass)
884     call chemmass1(ch4, ch4mass)
885    
886     ! === if hfc, pfc, and sf6 are included:
887     #ifdef INC_3GASES
888     ! === 032698
889     call chemmass1(hfc134a, hfc134amass)
890     call chemmass1(pfc, pfcmass)
891     call chemmass1(sf6, sf6mass)
892     ! ===
893     #endif
894    
895     !
896     ! === Calculating advection and eddy diffusion:
897     !
898     call chemadv0 (dt)
899    
900     !
901     ! === Calculate total n-s transport amount
902     ! === of cfc11 - temperary:
903     !
904     dth = 3600.0 ! for relexible setting of dt dt * 3.0
905     ! call chemtmp1 (dth,airmass0,p,pvv,cfc11)
906    
907     ! === Readjust mass of tracers, 1:
908     !
909     ! -----------------------------
910     ! Use tropospheric life time (yr) to calculate mass
911     ! loss of tracers and at the same time compensate
912     ! all the numerical loss back to tracer's mass
913     ! which equavalent to use adjcoe = 1.0 for
914     ! chemmass2.f
915     !
916     ! Chien Wang, September 12,1995
917     ! ------------------------------
918    
919     ! === 092595 update p
920    
921     call chemairmass(p4chem1)
922    
923     call chemmass6(46.0, 1.0, cfc11,cfc11mass)
924     call chemmass6(120.0,1.0, cfc12,cfc12mass)
925    
926     ! === 102596
927     ! === close tau type of ocean uptake of co2:
928    
929     call chemmass66(1.0, 1.0,zco2,zco2mass)
930    
931     call chemmass6(150.0,1.0,xn2o,xn2omass)
932     call chemmass2(1.0,ch4, ch4mass )
933    
934     ! === if hfc, pfc, and sf6 are included:
935     #ifdef INC_3GASES
936     ! === 032698
937     ! call chemmass6(14.6, 1.0,hfc134a, hfc134amass)
938     call chemmass2(1.0, hfc134a, hfc134amass)
939     call chemmass6(10000.0,1.0,pfc,pfcmass)
940     call chemmass6(3200.0, 1.0,sf6,sf6mass)
941     ! ===
942     #endif
943    
944     !
945     ! === Calculate tropospheric gaseous reactions
946     ! === every nhr_for_chem hours:
947     !
948     dt_chem_h = dth*float(nhr_for_chem)
949     if(ievenodd.eq.0) then
950     call chemtrop0(0, T, q, dt_chem_h, 1)
951    
952     c print *,'H2SO4 after chemtrop0 ',h2so4(1,33,1)
953     c print *,'SVIOD after chemtrop0 ',SVIOD(1,33,1)
954    
955     !
956     ! --- tmp output
957     !
958     cnt3hr(hrcnt) = cnt3hr(hrcnt) + 1
959     sfc3hro3(1:nlat,hrcnt) = sfc3hro3(1:nlat,hrcnt)
960     & + (o3(1,1:nlat,1))*29.0/48.0
961    
962     hrcnt = hrcnt + 1
963     if (hrcnt .gt. 8 ) hrcnt = 1
964     end if
965    
966     ievenodd = ievenodd + 1
967     if(ievenodd.ge.nhr_for_chem) ievenodd = 0
968    
969     !
970     ! === Calculating stratospheric processes:
971     ! === 092595
972     ! === adjust startospheric loss to whole global:
973     !
974     call chemmass1(cfc11, cfc11mass)
975     call chemmass1(cfc12, cfc12mass)
976     call chemmass1(xn2o, xn2omass)
977    
978     call chemstrat (dt)
979    
980     call chemmass2(1.0,cfc11,cfc11mass)
981     call chemmass2(1.0,cfc12,cfc12mass)
982     call chemmass2(1.0,xn2o, xn2omass)
983    
984     !
985     ! === Get total mass of chemically active species:
986     !
987     ! call chemairmass(p)
988     ! write(6,*)"P before 2nd", p
989     c print *,'H2SO4 after chemairmass ',h2so4(1,33,1)
990    
991     call chemmass1(cfc11, cfc11mass)
992     call chemmass1(cfc12, cfc12mass)
993     call chemmass1(xn2o, xn2omass)
994     call chemmass1(zco2, zco2mass)
995     call chemmass1(ch4, ch4mass)
996    
997     ! === if hfc, pfc, and sf6 are included:
998     #ifdef INC_3GASES
999     ! === 032698
1000     call chemmass1(hfc134a, hfc134amass)
1001     call chemmass1(pfc, pfcmass)
1002     call chemmass1(sf6, sf6mass)
1003     ! ===
1004     #endif
1005    
1006     #endif
1007    
1008     c print *,'H2SO4 before CONDSE ',h2so4(1,33,1)
1009     CALL CONDSE(mndriver) 141.
1010     c print *,'H2SO4 after CONDSE ',h2so4(1,33,1)
1011    
1012     #if ( defined CPL_CHEM )
1013     !
1014     ! === Calculate emission once per hour:
1015     !
1016     ! timeinhr=1./(365.*24.) !hourly emission
1017     ! call chememission(timeinhr)
1018     !
1019     !
1020     ! === Print hourly
1021     !
1022     ! call chemprt
1023     !
1024     ! ==============================================
1025     #endif
1026    
1027     if(HPRNT)then
1028     print *,' main after conds',' TAU=',TAU
1029     #include "PRNT.COM"
1030     endif
1031     #if ( !defined CLM )
1032     CALL PRECIP_LAND(mndriver) 142.
1033     if(HPRNT)then
1034     print *,' main after preci',' TAU=',TAU
1035     #include "PRNT.COM"
1036     endif
1037     #endif
1038     CALL CLOCKS (MNOW) 143.
1039     MINC=MLAST-MNOW 144.
1040     MCNDS=MCNDS+MLAST-MNOW 145.
1041     MLAST=MNOW 146.
1042     C 147.
1043     CALL DIAG9A (3) 148.
1044     C**** RADIATION, SOLAR AND THERMAL 149.
1045     if(HPRNT)then
1046     print *,' main before radia',' TAU=',TAU
1047     #include "PRNT.COM"
1048     print *,' IRAND=',IRAND
1049     endif
1050     #if ( defined PREDICTED_GASES || defined PREDICTED_AEROSOL )
1051     call radia_chem
1052     #else
1053     if(OBSFOR) then
1054     CALL RADIAGSO
1055     else
1056     CALL RADIA
1057     endif
1058     #endif
1059     if(HPRNT)then
1060     print *,' main after radia',' TAU=',TAU
1061     #include "PRNT.COM"
1062     endif
1063     CALL CLOCKS (MNOW) 151.
1064     MINC=MINC+MLAST-MNOW 152.
1065     MRAD=MRAD+MLAST-MNOW 153.
1066     MLAST=MNOW 154.
1067     if(HPRNT)then
1068     print *,' main before diag9a',' TAU=',TAU
1069     #include "PRNT.COM"
1070     endif
1071     C 155.
1072     CALL DIAG9A (4) 156.
1073     if(HPRNT)then
1074     print *,' main after diag9a',' TAU=',TAU
1075     #include "PRNT.COM"
1076     endif
1077     C**** SURFACE INTERACTION AND GROUND CALCULATION 157.
1078    
1079     #if ( defined CLM )
1080     if(HPRNT)then
1081     print *,' main before surf4clm',' TAU=',TAU
1082     #include "PRNT.COM"
1083     endif
1084     CALL SUR4CLM
1085     if(HPRNT)then
1086     print *,' main after surf4clm',' TAU=',TAU
1087     #include "PRNT.COM"
1088     endif
1089 jscott 1.4 i=1
1090 jscott 1.1 do j=1,jm
1091 jscott 1.4 pcpl4clm(i,j)=pcpl4clm(i,j)*prlnd2total(j,mndriver)
1092     pcpc4clm(i,j)=pcpc4clm(i,j)*prlnd2total(j,mndriver)
1093 jscott 1.1 enddo
1094     ! print *,' main after surf4clm',' TAU=',TAU
1095     ! print ('2(12f7.2,/,11f7.2,/)'),ps4clm,pcpl4clm,
1096     ! & pcpc4clm,tpr4clm,
1097     ! & tsl4clm,
1098     ! & qs4clm,ws4clm
1099     ! & ,us4clm,vs4clm,
1100     ! & dsw4clm,
1101     ! & dlw4clm,pco24clm
1102     ! & ,swinr4clm,swvis4clm
1103     #if ( defined DATA4TEM )
1104     c if(JYEAR.gt.20)then
1105     write (935),ps4clm,pcpl4clm,
1106     & pcpc4clm,tpr4clm,
1107     & tsl4clm,
1108     & qs4clm,ws4clm
1109     & ,us4clm,vs4clm,
1110     & dsw4clm,
1111     & dlw4clm,pco24clm
1112     & ,swinr4clm,swvis4clm
1113     c endif
1114     #endif
1115    
1116     ncallclm=ncallclm+1
1117     c print *,'before clm4mit2d ncallclm=',ncallclm
1118     if(HPRNT)then
1119     print *,' main before clm4mit2d',' TAU=',TAU
1120     #include "PRNT.COM"
1121     endif
1122     call clm4mit2d
1123     if(HPRNT)then
1124     print *,' main after clm4mit2d',' TAU=',TAU
1125     #include "PRNT.COM"
1126     endif
1127    
1128     c if(JYEAR.gt.20)then
1129 jscott 1.4 ! write (934),tau,tsoiclm,snwdclm,snwcclm,
1130     ! & lwuclm,tref2mclm,tflxclm,tgndclm,
1131     ! & lhfclm,shfclm,tauxclm,tauyclm,
1132     ! & asdirclm,aldirclm,asdifclm,aldifclm,
1133     ! & sroclm,ssrclm,glrclm
1134     ! &,h2olclm,h2oiclm
1135 jscott 1.1 c endif
1136     ! print *,' main after clm4mit2d',' TAU=',TAU
1137     ! print ('2(12f7.2,/,11f7.2,/)'),tsoiclm,snwdclm,snwcclm,
1138     ! & lwuclm,tref2mclm,tflxclm,tgndclm,
1139     ! & lhfclm,shfclm,tauxclm,tauyclm,
1140     ! & asdirclm,aldirclm,asdifclm,aldifclm
1141    
1142     CALL SURF_CLM
1143     CALL SURF_OCEAN
1144     CALL GR_CLM
1145     if(HPRNT)then
1146     print *,' main after surfc',' TAU=',TAU
1147     #include "PRNT.COM"
1148     endif
1149     #else
1150     ! CALL SURF_LAND
1151     ! CALL SURF_OCEAN
1152     CALL SURFCE ! 07/14/2006
1153     if(HPRNT)then
1154     print *,' main after surfc',' TAU=',TAU
1155     #include "PRNT.COM"
1156     endif
1157     CALL GRLAND 159.
1158     #endif
1159    
1160     #if ( defined OCEAN_3D || defined ML_2D)
1161     CALL GRFOROCEAN
1162     #else
1163     CALL GROCEAN(mndriver) 159.
1164     #endif
1165     c
1166     if(HPRNT)then
1167     print *,' main after groun',' TAU=',TAU
1168     #include "PRNT.COM"
1169     endif
1170     c print *,'H2SO4 before DRYCNV ',h2so4(1,33,1)
1171     CALL DRYCNV 160.
1172     c print *,'H2SO4 after DRYCNV ',h2so4(1,33,1)
1173     if(HPRNT)then
1174     print *,' main after drycn',' TAU=',TAU
1175     #include "PRNT.COM"
1176     endif
1177     CALL CLOCKS (MNOW) 161.
1178     MINC=MINC+MLAST-MNOW 162.
1179     MSURF=MSURF+MLAST-MNOW 163.
1180     MLAST=MNOW 164.
1181     CALL DIAG9A (5) 165.
1182     C**** STRATOSPHERIC MOMENTUM DRAG 166.
1183     CALL SDRAG(WLMMAX,JWLMMAX) 167.
1184     if(HPRNT)then
1185     print *,' main after sdrag',' TAU=',TAU
1186     #include "PRNT.COM"
1187     endif
1188     CALL CLOCKS (MNOW) 168.
1189     MINC=MINC+MLAST-MNOW 169.
1190     MSURF=MSURF+MLAST-MNOW 170.
1191     MLAST=MNOW 171.
1192     C 172.
1193     CALL DIAG9A (6) 173.
1194     MSRCE=MCNDS+MRAD+MSURF 174.
1195     PERCNT=100.*MSRCE/(MSTART-MNOW+1.E-5) 175.
1196     C**** SEA LEVEL PRESSURE FILTER 177.
1197     NFILTR= NDYN 177.5
1198     IF(MFILTR.LE.0.OR.MOD(NSTEP,NFILTR).NE.0) GO TO 500 178.
1199     IDACC(10)=IDACC(10)+1 179.
1200     if(HPRNT)then
1201     print *,' main before filtr',' TAU=',TAU
1202     #include "PRNT.COM"
1203     endif
1204     C
1205     C 180.
1206     C ****************
1207     if(LHORDIF)then
1208     DTDIF=3600.
1209     if(JM.eq.24)then
1210     CALL HORDIFF(DTDIF)
1211     else if(JM.eq.46)then
1212     CALL HORDIFFALL(DTDIF)
1213     else
1214     print *,' Wromg JM=',JM
1215     stop
1216     endif
1217     end if
1218     C ****************
1219     c CALL FILTER 181.
1220     C
1221     #if ( defined CPL_CHEM )
1222     !
1223     ! === Readjust total mass of tracers:
1224     !
1225     call chemairmass(p)
1226    
1227     call chemmass2(1.00,cfc11,cfc11mass)
1228     call chemmass2(1.00,cfc12,cfc12mass)
1229     call chemmass2(1.00,xn2o ,xn2omass )
1230     call chemmass2(1.00,zco2,zco2mass)
1231     call chemmass2(1.00,ch4, ch4mass )
1232    
1233     ! === if hfc, pfc, and sf6 are included:
1234     #ifdef INC_3GASES
1235     ! === 032698
1236     call chemmass2(1.0,hfc134a, hfc134amass)
1237     call chemmass2(1.0,pfc, pfcmass)
1238     call chemmass2(1.0,sf6, sf6mass)
1239     ! ===
1240     #endif
1241    
1242     !
1243     ! === Accumulative calculation prepared for
1244     ! === carrying out monthly average:
1245     !
1246     call chemmonth1
1247     !
1248     #endif
1249    
1250     CALL CLOCKS (MNOW) 182.
1251     MDYN=MDYN+MLAST-MNOW 183.
1252     MLAST=MNOW 184.
1253     C 185.
1254     CALL DIAG9A (7) 186.
1255     C**** 187.
1256     C**** UPDATE MODEL TIME AND CALL DAILY IF REQUIRED 188.
1257     C**** 189.
1258     500 NSTEP=NSTEP+NDYN 190.
1259     ITAU=(NSTEP+NSTEP0)*IDTHR 191.
1260 jscott 1.7 cJRS fix to DFLOAT 8/2/07
1261     TAU=DFLOAT(ITAU)/XINT 192.
1262 jscott 1.1 IDAY=1+ITAU/I24 193.
1263     TOFDAYPR=TOFDAY+1.00
1264     TOFDAY=(ITAU-(IDAY-1)*I24)/XINT 194.
1265     IF(.NOT.EVENT(24.)) GO TO 550 195.
1266     C 196.
1267     do J=1,JM0
1268     TSURFW(J)=TSURFD(J)
1269 jscott 1.4 TLANDW(J)=TLANDD(J)
1270 jscott 1.1 enddo
1271    
1272     JDATECLM=JDATE
1273     JDATEPR=JDATE
1274     JMONTHPR=JMONTH
1275     JYEARPR=JYEAR
1276     CALL DAILY_NEW 197.
1277     c print *,' AFTER DAILY_NEW IDAY=',IDAY,' IYEAR=',IYEAR
1278     c print *,' JYEAR=',JYEAR,' JDAY=',JDAY
1279     #if( !defined OCEAN_3D && !defined ML_2D )
1280     CALL DAILY_OCEAN
1281     c print *,' AFTER DAILY_OCEAN IDAY=',IDAY,' IYEAR=',IYEAR
1282     c print *,' JYEAR=',JYEAR,' JDAY=',JDAY
1283     c print *,' JDATE=',JDATE,' JMONTH=',JMONTH
1284     #endif
1285     if(JDATE.eq.100)then
1286     print *,JDATE,JMONTH,JYEAR
1287     print *,' main after daily'
1288     print *,' T1 ocean'
1289     print 5001,(ODATA(1,j,1),j=1,JM0)
1290     print *,' T2 ocean'
1291     print 5001,(ODATA(1,j,4),j=1,JM0)
1292     print *,' T3 ocean'
1293     print 5001,(ODATA(1,j,5),j=1,JM0)
1294     print *,' sea ice'
1295     print 5002,(ODATA(1,j,2),j=1,JM0)
1296     endif
1297     CALL CLOCKS (MNOW) 198.
1298     MELSE=MELSE+(MLAST-MNOW) 199.
1299     MLAST=MNOW 200.
1300     NDAILY=SDAY/DT 201.
1301    
1302     #if ( defined CPL_CHEM )
1303     !
1304     ! === Calculate air mass, second step:
1305     !
1306     i=1
1307     call chemairmass(p)
1308    
1309     !
1310     ! === Calculate emission once per day:
1311     !
1312     c zxy=0.
1313     c do j=1,jm
1314     c zxy=zxy+zco2(1,j,1)
1315     c & *28.97296245/44.0*1.e-3
1316     c enddo
1317     c print *,' CO2 before emission ',zxy/jm
1318    
1319     timeinday=1./(365.) !daily emission
1320     call chememission(timeinday)
1321    
1322     c zxy=0.
1323     c do j=1,jm
1324     c zxy=zxy+zco2(1,j,1)
1325     c & *28.97296245/44.0*1.e-3
1326     c enddo
1327     c print *,' CO2 after emission ',zxy/jm
1328     !
1329     #endif
1330    
1331     C 202.
1332     CALL DIAG9A (8) 203.
1333     #if ( !defined OCEAN_3D && !defined ML_2D )
1334     IF(KOCEAN.EQ.1) THEN
1335     DO 540 J=1,JM 203.11
1336     DO 540 I=1,IM 203.12
1337     AIJ(I,J,59)=AIJ(I,J,59)+ODATA(I,J,4) 203.13
1338     540 AIJ(I,J,60)=AIJ(I,J,60)+ODATA(I,J,5) 203.14
1339     if(TRANSR)then
1340     if(LMO.eq.11) then
1341     CALL ODIFS
1342     elseif(LMO.eq.12) then
1343     CALL ODIFS12
1344     else
1345     Print *,' Wromng LMO',LMO
1346     stop
1347     endif
1348     NODIFS=NODIFS+1
1349     do 5283 j=1,JM0
1350     do 5283 k=1,lmo
1351     RTGOAV(j,k)=RTGOAV(j,k)+RTGO(1,j,k)
1352     5283 continue
1353     endif
1354     C**** RESTRUCTURE THE OCEAN LAYERS AND ELIMINATE SMALL ICE BERGS 203.15
1355     CALL OSTRUC 203.16
1356     if(JDATE.eq.41)then
1357     print *,JDATE,JMONTH,JYEAR
1358     print *,' main after ostruc'
1359     print *,'Z1O'
1360     print 5001,(Z1O(1,j),j=1,JM0)
1361     print *,' T1 ocean'
1362     print 5001,(ODATA(1,j,1),j=1,JM0)
1363     print *,' T2 top ocean'
1364     print 5001,((2.*ODATA(1,j,4)-ODATA(1,j,5)),j=1,JM0)
1365     print *,' T2 ocean'
1366     print 5001,(ODATA(1,j,4),j=1,JM0)
1367     print *,' T3 ocean'
1368     print 5001,(ODATA(1,j,5),j=1,JM0)
1369     endif
1370     ENDIF ! KOCEAN
1371     #endif
1372     #if ( defined ML_2D )
1373     if(TRANSR)then
1374     NODIFS=NODIFS+1
1375     do 5283 j=1,JM0
1376     do 5283 k=1,lmo
1377     RTGOAV(j,k)=RTGOAV(j,k)+RTGO(1,j,k)
1378     5283 continue
1379     endif
1380     #endif
1381    
1382     #if ( defined CPL_OCEANCO2 )
1383    
1384     C For OCM or 3D ocean model with carbon cycle
1385     #ifdef PREDICTED_GASES
1386     c -------
1387     c 102596
1388     c air co2 mixing ratio goes to ocean carbon model:
1389     c
1390     do j=1,jm
1391     pC_atm(j)=zco2(1,j,1)
1392     & *28.97296245/44.0*1.e-9
1393     !ppb(m) to kg per volume base
1394 jscott 1.4
1395     atm_co2(j)=pC_atm(j)*1.e6
1396    
1397 jscott 1.1 enddo ! j
1398     c
1399     c -------
1400     #else
1401    
1402     do j=1,jm
1403     pC_atm(j)=atm_co2(j)*1.e-6
1404     enddo ! j
1405    
1406     #endif
1407    
1408     #if ( defined ML_2D )
1409     CB Gary CO2 uptake by the ocean
1410    
1411     do j=1,jm
1412     if(NWMGEA(J).gt.0)then
1413     WSAV(J)=WSAV(J)/NWMGEA(J)
1414     NWMGEA(J)=0.
1415     else
1416     WSAV(J)=0.
1417     end if
1418     end do ! j
1419    
1420     do j=1,jm
1421     tggary(j)=ODATA(1,j,1)+273.16
1422     wsgary(j)=WSAV(J)
1423     WSAV(j)=0.
1424     arsrf(j)=areaml(j)*(1.-ODATA(1,j,2))
1425     DEPTHML(j)=ZOAV(J)
1426     co24ocnan(j)=co24ocnan(j)+pC_atm(j)*1.e6
1427     enddo ! j
1428     c print *,'CO2 for 2D ocean'
1429     c print *,JYEAR,JMONTH
1430     c print *,'co2=',pC_atm(27)*1.e6,' ws=',wsgary(27)
1431     c print *,'tem=',tggary(27)
1432     c print '12f7.1,/,2(11f7.1,/,),12f7.1',(pC_atm(j)*1.e6,j=1,jm)
1433     c print '12f7.1,/,2(11f7.1,/,),12f7.1',(rco2(j,1),j=1,jm)
1434     c ncallgary=ncallgary+1
1435 jscott 1.4 ! 10/28/06
1436     ! call carb_mxdlyr_chem(focean)
1437     ! call carb_airsea_flx
1438     !
1439     ! 3D ocean chemistry
1440     call carb_chem_ocmip(focean)
1441     call carb_airsea_flx(dtco2)
1442     ! 3D ocean chemistry
1443     ! 10/28/06
1444 jscott 1.1 c print *,'FCO2 ncallgary=',ncallgary
1445     c print '12f7.1,/,2(11f7.1,/,),12f7.1',
1446     c & (fluxco2(j)*12.e-15*365.,j=1,jm)
1447     #endif
1448    
1449     C For ocean carbon model
1450     c Annual oceanic CO2 uptake
1451     do j=1,jm
1452     OCUPT=OCUPT+fluxco2(j)
1453     enddo
1454 jscott 1.5 ! print *,' OCUPT=',OCUPT*12.e-15
1455 jscott 1.1
1456     #if ( defined CPL_CHEM )
1457     !
1458     ! === Calculate ocean uptake of CO2
1459     ! === once per day:
1460     !
1461     i=1
1462     call chemairmass(p) !update airmass
1463    
1464     call chemoceanco2(fluxco2)
1465     !
1466     #endif
1467    
1468 jscott 1.4 do j=1,jm
1469     fluxco2(j)=0.0
1470     enddo
1471 jscott 1.1 #endif
1472    
1473     #if ( defined CPL_TEM )
1474     !#if ( defined CLM )
1475     c print *,'JDATE for TEM=',JDATECLM
1476     do j=1,jm
1477     if(npred4tem(j).gt.0)then
1478     c pred4tem(j)=pred4tem(j)/npred4tem(j)
1479     ewvd4tem(j)=ewvd4tem(j)/npred4tem(j)
1480     pre4tem(J)=pre4tem(J)+pred4tem(j)
1481     endif
1482     c
1483     if(nt2md4tem(j).gt.0)then
1484     t2md4tem(j)=t2md4tem(j)/nt2md4tem(j)
1485     endif
1486     temp4tem(j)=temp4tem(j)+t2md4tem(j)
1487     dtem4tem(JDATECLM,j)=t2md4tem(j)
1488     c
1489     if(nradd4tem(j).gt.0)then
1490     cldd4tem(j)=cldd4tem(j)/ncldd4tem(j)
1491     swtd4tem(j)=swtd4tem(j)/nradd4tem(j)
1492     swsd4tem(j)=swsd4tem(j)/nradd4tem(j)
1493     endif
1494     sws4tem(j)=sws4tem(j)+swsd4tem(j)
1495     c
1496     enddo
1497     c
1498     do j=1,jm
1499     pred4tem(j)=0.0
1500     ewvd4tem(j)=0.0
1501     t2md4tem(j)=0.0
1502     cldd4tem(j)=0.0
1503     swtd4tem(j)=0.0
1504     swsd4tem(j)=0.0
1505     npred4tem(j)=0
1506     ncldd4tem(j)=0
1507     nradd4tem(j)=0
1508     nt2md4tem(j)=0.
1509     enddo
1510     #endif
1511    
1512     #if ( defined CPL_OCEANCO2 && defined ML_2D )
1513     C For OCM
1514    
1515 jscott 1.4 ! dtco2=3600.*24.
1516     ! cfkvct=1.0
1517     ! if (JYEAR.ge.1991) then begin
1518     ! if (JYEAR.le.2100) then begin
1519     ! cfkvct=(1.0*(2100-JYEAR)+0.25*(JYEAR-1990))/110.
1520     ! esle
1521     ! cfkvct=0.25
1522     ! endif
1523     ! endif
1524     ! do j=1,jm
1525     ! edzcatr(j)=cfkvct*edzcar(j)
1526     ! enddo
1527     call diffusco2(lmo,jm,dtco2,0.5,edzcart,depthml,focean,
1528 jscott 1.1 & dzg,dzog,rco2)
1529     call hdocean(rco2,focean,dxv,dyv,DXYP,depthml,edohd,dtco2)
1530     call avegary
1531     CB Gary CO2 uptake by the ocean
1532     #endif
1533    
1534     #if ( defined CPL_CHEM) && ( defined CPL_TEM )
1535     C take into accout land uptake form TEM for previous month
1536     do j=1,jm
1537     fluxnep(j)=aduptd+1.e-3*temco2(j)/NDAYMN(mndriver)
1538     c Annual TEM CO2 uptake
1539     TEMUPTANN=TEMUPTANN+fluxnep(j)
1540     enddo
1541     if(jdate.eq.1)then
1542     print *,'Monthly TEM uptake'
1543     c print *,mndriver,adupt,temuptann-temup0
1544     temup0=temuptann
1545     endif
1546     C
1547     c
1548     i=1
1549     call chemairmass(p) !update airmass
1550    
1551     call chemtemco2(fluxnep)
1552     C
1553     c
1554     !
1555     #endif
1556    
1557     c End of month
1558     if(JDATE.eq.1)then
1559    
1560     #if ( defined CPL_CHEM )
1561     !
1562     ! === Calculating monthly averaged mixing ratios:
1563     !
1564     call chemmonth2
1565     c print *,'Atmosphere after chemmonth2'
1566    
1567     ! === Calculate and print monthly n-s transport
1568     ! === of cfc11:
1569     !
1570     ! call chemtmp2
1571     !
1572     do nhr = 1,8
1573     sfc3hro3(1:nlat,nhr) = sfc3hro3(1:nlat,nhr)
1574     & /float(cnt3hr(nhr))
1575     #if ( defined CPL_TEM )
1576     o34tem(nhr,1:nlat)=sfc3hro3(1:nlat,nhr)
1577     #endif
1578     end do
1579     cwrite(124)sfc3hro3
1580     cnt3hr(1:8) = 0
1581     sfc3hro3(1:nlat,1:8) = 0.0
1582    
1583     ! === Writing rawdata every month:
1584     !
1585     ! call chemprt !closed 032697
1586    
1587     call chembudget (p)
1588     print *,' Atmosphre after chembudget mymonth=',mymonth
1589     ! === 09/26/94
1590     ! === Reset year and month index:
1591     !
1592     mymonth = mymonth + 1
1593     if(mymonth.gt.12)then
1594     myyear = myyear +1
1595     ! myyear = min(myyear,nchemyr)
1596     myyear = min(myyear,myyearlast)
1597     mymonth = 1
1598     ! endif ! 27/8/2005
1599    
1600     ! === 092295
1601     ! === write rawdata for possible renew run
1602     ! === at end of each month:
1603     ! === at end of each year: 27/8/2005
1604     !
1605     rewind 178
1606     print *,'For chem restart ',myyear,mymonth
1607     write(178)myyear,mymonth,airmass,
1608     & cfc11,cfc110,
1609     & cfc11m,
1610     & cfc11sd,
1611     & cfc12,cfc12m,
1612     & cfc12sd,
1613     & xn2o ,xn2om ,
1614     & xn2osd,
1615     & hfc134a,hfc134am,
1616     & pfc ,pfcm,
1617     & sf6 ,sf6m,
1618     & bcarbon,bcm,
1619     & ocarbon,ocm,
1620     & atomo ,
1621     & o1d ,
1622     & o3 ,o3m ,
1623     & co ,com ,
1624     & zco2 ,zco2m,
1625     & atomh ,
1626     & ho ,
1627     & ho2 ,hoxm ,
1628     & h2o2 ,
1629     & xno ,
1630     & xno2 ,xnoxm,
1631     & xno3 ,
1632     & xn2o5 ,xnoym,
1633     & hno3 ,
1634     & ch4 ,ch4m ,
1635     & ch3 ,
1636     & cho ,
1637     & ch2o ,
1638     & ch3o ,
1639     & ch3o2 ,
1640     & ch3o2h,
1641     & so2 ,so2m ,
1642     & hoso2 ,
1643     & so3 ,
1644     & h2so4 ,h2so4m,
1645     & sviod ,sviodm
1646     rewind 178
1647    
1648     endif
1649     #endif
1650    
1651     #if ( defined CPL_TEM )
1652     !#if ( defined CLM )
1653     do J=1,JM
1654     c
1655     #ifdef PREDICTED_GASES
1656     co24tem(j)=zco2(1,j,1)
1657     & *28.97296245/44.0*1.e-3
1658     !ppm(m) to kg per volume base
1659     #else
1660     co24tem(j)=atm_co2(j)
1661    
1662     #endif
1663     c
1664     enddo
1665     #endif
1666    
1667     endif ! end of month
1668     550 continue
1669     C END of DAilY
1670    
1671    
1672     C CALL CHECKT (11) 203.17
1673     CALL CLOCKS (MNOW) 203.18
1674     MSURF=MSURF+(MLAST-MNOW) 203.19
1675     MLAST=MNOW 203.2
1676     C**** 204.
1677     C**** WRITE INFORMATION ONTO A TAPE EVERY USET HOURS 205.
1678     C**** 206.
1679     IF(USET.LE.0.) GO TO 600 207.
1680     IF(.NOT.EVENT(USET)) GO TO 600 208.
1681     C COMPUTATIONS FOR XXXXXX 209.
1682     WRITE (520) TAU,XXXXXX 210.
1683     CALL CLOCKS (MNOW) 211.
1684     MINC=MLAST-MNOW 212.
1685     MELSE=MELSE+MINC 213.
1686     PERCNT=100.*MELSE/(MSTART-MNOW+1.E-5) 214.
1687     c WRITE (6,910) MINC,MELSE,PERCNT,TAU 215.
1688     C**** 216.
1689     C**** CALL DIAGNOSTIC ROUTINES 217.
1690     C**** 218.
1691     600 IF(MOD(NSTEP-NDYN,NDA4).EQ.0) CALL DIAG4A 219.
1692     C 220.
1693     IF(NDPRNT(1).GE.0) GO TO 610 221.
1694     C**** PRINT CURRENT DIAGNOSTICS (INCLUDING THE INITIAL CONDITIONS) 222.
1695     IF(KDIAG(1).LT.9) CALL DIAG1(NOCLM) 223.
1696     IF(KDIAG(2).LT.9) CALL DIAG2 224.
1697     IF(KDIAG(7).LT.9) CALL DIAG7P 225.
1698     IF(KDIAG(3).LT.9) CALL DIAG3 226.
1699     C 227.
1700     C 228.
1701     IF(KDIAG(4).LT.9) CALL DIAG4 229.
1702     IF(TAU.LE.TAUI+DTHR*(NDYN+.5)) CALL DIAGKN 230.
1703     NDPRNT(1)=NDPRNT(1)+1 231.
1704     C 690 changed to 691 02/21/2003
1705     610 IF(.NOT.EVENT(24.)) GO TO 691 232.
1706     C**** PRINT DIAGNOSTIC TIME AVERAGED QUANTITIES ON NDPRNT-TH DAY OF RUN 233.
1707     DO 620 K=1,13 234.
1708     IF (JDAY.EQ.NDPRNT(K)) GO TO 630 235.
1709     620 CONTINUE 236.
1710     GO TO 640 237.
1711     c 630 WRITE (6,920) 238.
1712     630 continue
1713     IF(KDIAG(1).LT.9) CALL DIAG1(NOCLM) 239.
1714     IF(KDIAG(2).LT.9) CALL DIAG2 240.
1715     IF(KDIAG(7).LT.9) CALL DIAG7P 241.
1716     IF(KDIAG(3).LT.9) CALL DIAG3 242.
1717     C 243.
1718     C 244.
1719     C IF(KDIAG(6).LT.9) CALL DIAG6 245.
1720     IF(KDIAG(4).LT.9) CALL DIAG4 246.
1721     IF(KDIAG(8).LT.9) CALL DIAG8 (IPFLAG) 247.
1722     C**** THINGS TO DO BEFORE ZEROING OUT THE ACCUMULATING ARRAYS 248.
1723     C**** (NORMALLY DONE AT THE END OF A MONTH) 249.
1724     640 DO 650 K=1,13 250.
1725     IF(JDAY.EQ.NDZERO(K)) GO TO 660 251.
1726     650 CONTINUE 252.
1727     GO TO 690 253.
1728     C**** PRINT THE KEY DIAGNOSTICS 254.
1729     660 CONTINUE 255.
1730     CALL DIAGKN 255.
1731     C**** PRINT AND ZERO OUT THE TIMING NUMBERS 256.
1732     CALL CLOCKS (MNOW) 257.
1733     MDIAG=MDIAG+(MLAST-MNOW) 258.
1734     MLAST=MNOW 259.
1735     TOTALT=.01*(MSTART-MNOW) 260.
1736     PDYN=MDYN/TOTALT 261.
1737     PCDNS=MCNDS/TOTALT 262.
1738     PRAD=MRAD/TOTALT 263.
1739     PSURF=MSURF/TOTALT 264.
1740     PDIAG=MDIAG/TOTALT 265.
1741     PELSE=MELSE/TOTALT 266.
1742     DTIME=24.*TOTALT/(60.*(TAU-TAU0)) 267.
1743     c WRITE (6,909) DTIME,PDYN,PCDNS,PRAD,PSURF,PDIAG,PELSE 268.
1744     MDYN=0 269.
1745     MCNDS=0 270.
1746     MRAD=0 271.
1747     MSURF=0 272.
1748     MDIAG=0 273.
1749     MELSE=0 274.
1750     MSTART=MNOW 275.
1751     if(TRANSR)then
1752     do 5284 j=1,JM0
1753     do 5284 k=1,lmo
1754     RTGOAV(j,k)=RTGOAV(j,k)/NODIFS
1755     5284 continue
1756     c print *,'ATM RTGOAV monthly'
1757     c print *,(RTGOAV(J,1),j=1,jm)
1758     endif ! TRANSR
1759     SPGAV=srps(jm+3)/nsrps+PTOP
1760     do 5285 j=1,JM+3
1761     GBUDG(j,38,1)=(srps(j)/nsrps+PTOP)*1013./SPGAV
1762     5285 continue
1763     c do 5287 j=1,JM+3
1764     c GBUDG(j,38,1)=GBUDG(j,37,1)*1013./GBUDG(jm+3,37,1)
1765     c5287 continue
1766 jscott 1.4 ! print *,'FRMDICE'
1767     ! print '6(1PE12.4)',FRMDICE
1768 jscott 1.1 ENKE=0.0
1769     ENPT=0.0
1770     do ii=1,4
1771     ENKE=ENKE+SPECA(1,19,ii)
1772     ENPT=ENPT+SPECA(1,20,ii)
1773     enddo
1774     c print *,'ENKE=',ENKE,' ENTP=',ENPT,' ENTT=',ENKE+ENPT
1775     print *,'ENKE=',QTABLE(JMP3,19,10),' ENTP=',QTABLE(JMP3,20,10),
1776     & ' ENTT=',QTABLE(JMP3,19,10)+QTABLE(JMP3,20,10)
1777     print *,'WLMMAX=',WLMMAX,' JWLMMAX=',JWLMMAX
1778     c print *,'AJ(*,37)'
1779     c print *,(AJ(J,37),j=1,jm)
1780     c print *,'AJ(*,28)'
1781     c print *,(AJ(J,28),j=1,jm)
1782     c IF(USEP.LE.0.) GO TO 680 276.
1783     C**** WRITE SELECTED DIAGNOSTICS ONTO A DISK DATA SET FOR PLOTTING 277.
1784     c IF(TAU.LE.TAUI+1080.) GO TO 675 278.
1785     c 670 READ (16) TAUX 279.
1786     c IF(TAU.GT.TAUX+1080.) GO TO 670 280.
1787     675 WRITE (546) AEXP4,JDATE,JMONTH,JYEAR,JDATE0,JMNTH0,JYEAR0, 281.
1788     * GBUDG,QMAPS,QTABLE,INQTAB,J1QT,INQMAP,RTGOAV
1789     print *,'From atm write(546) ',JMNTH0,JYEAR0,' ',JMONTH,JYEAR
1790     c print *,0.1*GBUDG(1,26,2),0.1*GBUDG(1,35,2)
1791     nwr=nwr+1
1792     690 continue
1793     if(JDAY.eq.1)then
1794     rewind 547
1795     write(547)AEXP,nwr
1796     rewind 547
1797     print *,'From atm write(547) ',AEXP,nwr
1798     endif
1799     if(wr25.and.JDAY.eq.1)then
1800     c Write a restart file once a year
1801     print *,'Write a restart file once a year.'
1802    
1803     #if ( defined CPL_OCEANCO2 && defined ML_2D )
1804     C Data for possible restart for OCM
1805     write(369)jyear-1,vdfocm
1806     write(369)Hg
1807     write(369)Rco2
1808     rewind 369
1809     #endif
1810    
1811     print *,' KDISK0=',KDISK0
1812     CALL RFINAL (IRAND)
1813     REWIND KDISK0
1814     if(TRANSR)then
1815     WRITE(KDISK0) AEXP,TAU,JC,C,RC,KEYNR,U,V,T,P,Q,ODATA,GDATA,
1816     * BLDATA,
1817     * RQT,SRHR,TRHR,(AJ(K,1),K=1,KACC),TAU,TSSFC,CKS,CKN,WMGE,TPRIM2
1818 jscott 1.4 * ,MRCHT,TRSURF,SRSURF,TLANDW,TSURFW,DWAV0
1819 jscott 1.1 * ,TG3M,RTGO,STG3,DTG3
1820     else
1821     WRITE(KDISK0) AEXP,TAU,JC,C,RC,KEYNR,U,V,T,P,Q,ODATA,GDATA,
1822     * BLDATA,
1823     * RQT,SRHR,TRHR,(AJ(K,1),K=1,KACC),TAU,TSSFC,CKS,CKN,WMGE,TPRIM2
1824 jscott 1.4 * ,MRCHT,TRSURF,SRSURF,TLANDW,TSURFW,DWAV0
1825 jscott 1.1 c print *,' TSURFT'
1826     c print 5001,TSURFT
1827     c print *,' TSURFW'
1828     c print 5001,TSURFW
1829     endif
1830     REWIND KDISK0
1831     KDISK=3-KDISK
1832     KDISK0=500+KDISK
1833     end if
1834     C 690 changed to 691 02/21/2003
1835     680 IF(KCOPY.LE.0) GO TO 691 284.
1836     C**** WRITE A COPY OF THE FINAL RESTART DATA SET ONTO DISK 285.
1837     CALL RFINAL (IRAND) 285.5
1838     print *,' after 680'
1839     print *,' TAU=',TAU,' IRAND=',IRAND
1840     IF(KCOPY.GT.99) GO TO 687 286.
1841     683 READ (KCOPY) TAUX 286.5
1842     IF(TAU.GT.TAUX+3240.) GO TO 683 287.
1843     685 WRITE (KCOPY) TAU,JC,C,RC,KEYNR,U,V,T,P,Q,ODATA,GDATA,BLDATA, 287.5
1844     * RQT,SRHR,TRHR,(AJ(K,1),K=1,KACC),TAU,TSSFC,CKS,CKN 288.
1845 jscott 1.4 * ,WMGE,TPRIM2,MRCHT,TRSURF,SRSURF,TLANDW,TSURFW,DWAV0
1846 jscott 1.1 * ,TG3M,RTGO,STG3,DTG3
1847     REWIND KCOPY 288.5
1848     C 690 changed to 691 02/21/2003
1849     GO TO 691 289.
1850     687 KCOPY=KCOPY-100 289.5
1851     GO TO 685 289.6
1852     C**** TIME FOR CALLING DIAGNOSTICS 290.
1853     C 690 changed to 691 02/21/2003
1854     691 CALL CLOCKS (MNOW) 291.
1855     MDIAG=MDIAG+(MLAST-MNOW) 292.
1856     MLAST=MNOW 293.
1857     780 IF(TAU.LE.TAUI+DTHR*(NDYN+.5)) GO TO 120 294.
1858     GO TO 100 295.
1859     C**** 296.
1860     C**** END OF MAIN LOOP 297.
1861     C**** 298.
1862     C**** RUN TERMINATED BECAUSE SENSE SWITCH 6 WAS TURNED ON 299.
1863     800 WRITE (6,904) 300.
1864     IF(EVENT(TAUT)) GO TO 820 301.
1865     CALL RFINAL (IRAND) 302.
1866     print *,' after 800'
1867     print *,' TAU=',TAU,' IRAND=',IRAND
1868     if(wr25) then
1869     REWIND KDISK0 303.
1870     WRITE(KDISK0) AEXP,TAU,JC,C,RC,KEYNR,U,V,T,P,Q,ODATA,GDATA,BLDATA, 304.
1871     * RQT,SRHR,TRHR,(AJ(K,1),K=1,KACC),TAU,TSSFC,CKS,CKN,WMGE,TPRIM2 305.
1872 jscott 1.4 * ,MRCHT,TRSURF,SRSURF,TLANDW,TSURFW,DWAV0
1873 jscott 1.1 end if
1874     C WRITE (6,908) 306.
1875     C**** RUN TERMINATED BECAUSE IT REACHED TAUE (OR SS6 WAS TURNED ON) 307.
1876     c 820 WRITE (6,905) TAU,IDAY,TOFDAY 308.
1877     820 continue
1878    
1879     #if ( defined OCEAN_3D || defined ML_2D )
1880     C DTATM time step of atm model in hours
1881     C precip and evap in mm/day or kg/m**2/day
1882     do j=1,jm0
1883 jscott 1.8 Cjrs #if ( defined OCEAN_3D && defined CPL_OCEANCO2 )
1884     #ifdef OCEAN_3D
1885     !jrs ncallatm=ncallatm+1
1886 jscott 1.4 ! 020107
1887     ! co24ocean(j)=pC_atm(j)*1.e6
1888 jscott 1.8 ! jrs give CO2 even if ocn carbon off
1889 jscott 1.4 co24ocean(j)=atm_co2(j)
1890 jscott 1.8 # ifdef CPL_OCEANCO2
1891 jscott 1.1 co24ocnan(j)=co24ocnan(j)+co24ocean(j)
1892 jscott 1.8 # endif
1893 jscott 1.1 #endif
1894     #ifdef ML_2D
1895 jscott 1.4 cjrs block only MD_2D
1896 jscott 1.1 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 jscott 1.4 Cjrs surfpr(j)=surfpr(j)/(DTATM/24.)
1908 jscott 1.1 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 jscott 1.4 cjrs change runoff to new name arunoff
1946 jscott 1.1 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 jscott 1.4 c jrs runoff->arunoff
1976 jscott 1.1 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 jscott 1.4 c jrs only ML_2D
1990 jscott 1.1 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 jscott 1.4 c jrs only ML_2D
2039 jscott 1.1 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 jscott 1.4 Cjrs change this block to only ML_2D
2062 jscott 1.1 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 jscott 1.4 cjrs print *,'ncallclm=',ncallclm
2078 jscott 1.1 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