/[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.4 - (hide annotations) (download)
Mon Apr 23 21:20:17 2007 UTC (18 years, 3 months ago) by jscott
Branch: MAIN
Changes since 1.3: +122 -71 lines
bring igsm atmos code up to date

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