/[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.12 - (hide annotations) (download)
Thu Sep 17 15:48:38 2009 UTC (15 years, 10 months ago) by jscott
Branch: MAIN
CVS Tags: HEAD
Changes since 1.11: +7 -17 lines
new routine for reading in eppa emissions

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

  ViewVC Help
Powered by ViewVC 1.1.22