/[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.11 - (hide annotations) (download)
Tue Sep 1 21:56:40 2009 UTC (15 years, 10 months ago) by jscott
Branch: MAIN
Changes since 1.10: +6 -2 lines
change atmos lifetime of CO2

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

  ViewVC Help
Powered by ViewVC 1.1.22