--- MITgcm_contrib/jscott/igsm/src/atmosphere.F 2006/08/22 20:25:52 1.2 +++ MITgcm_contrib/jscott/igsm/src/atmosphere.F 2009/09/01 21:56:40 1.11 @@ -1,3 +1,5 @@ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/jscott/igsm/src/atmosphere.F,v 1.11 2009/09/01 21:56:40 jscott Exp $ +C $Name: $ #include "ctrparam.h" @@ -108,40 +110,42 @@ logical odifcarbon #if ( defined CLM ) -#include "CLM.COM" -!#include "TEM.COM" -#if ( defined CPL_TEM ) +# include "CLM.h" +# if ( defined CPL_TEM ) C For TEM -#include "TEM.COM" -#endif +# include "TEM.h" +# endif #endif + #if ( defined CPL_OCEANCO2 ) -#include "OCM.COM" +#include "OCM.h" common /Garyflux/pC_atm(jm0),wind_amp,fluxco2(jm0) # if ( defined ML_2D ) common/Garyclim/tggary(jm0),wsgary(jm0),areaml(jm0),arsrf(jm0) - common/Garydiff/depthml(jm0),edzon(jm0),dzg(lmo),dzog(lmo-1), + common/Garydiff/depthml(jm0),edzcar(jm0),dzg(lmo),dzog(lmo-1), &Rco2(jm0,lmo),edohd(lmo),zg(lmo),focean(jm0) common /Garychem/Hg(jm0) common/qfl/QFLUX(JM0,0:13),ZOAV(JM0),QFLUXT(JM0) common /Garyvdif/iyearocm,vdfocm,acvdfc - common /Garyvlog/odifcarbon + common /Garyvlog/odifcarbon,ocarcont + common /Garykvct/cfkvct,edzcart(jm0) # endif #endif - INTEGER dtatm, mndriver !routine arguments + INTEGER dtatm, mndriver !routine arguments jrs #if ( defined OCEAN_3D || defined ML_2D ) #include "AGRID.h" -C#include "HRD4OCN.COM" +Cjrs elimated COM file/moved elsewhere#include "HRD4OCN.COM" dimension oimeltt(jm0),dhdtav(jm0),devdtav(jm0) #endif C **** CLEAR SKY common/clrsk/CLEAR(JM0),NCLR(JM0),AJCLR(JM0,12),BJCLR(JM0,12), * CJCLR(JM0,12) - common/TSUR/TSURFC(JM0,0:13),TSURFT(JM0),TSURFD(JM0),DTSURF(JM0) - real TSURFW(JM0) +! common/TSUR/TSURFC(JM0,0:13),TSURFT(JM0),TSURFD(JM0),DTSURF(JM0) +#include "TSRF.COM" + real TSURFW(JM0),TLANDW(JM0) integer CLEAR common /ATCO2/atm_co2(jm0) @@ -205,11 +209,6 @@ IPFLAG=0 15. C CALL ENQJOB 16. CALL INPUT 17. - print *,"After input" - print *,"TSURFD" - print *,TSURFD - print *,"TSURFT" - print *,TSURFT #if ( defined CPL_CHEM ) ! @@ -257,6 +256,7 @@ #if ( defined CPL_OCEANCO2 && defined ML_2D ) odifcarbon=.true. wind_amp=1. + dtco2=3600.*24. c ncallgary=0 do j=1,jm areaml(j)=dxyp(j)*(1-FDATA(1,J,2)) @@ -264,7 +264,8 @@ DEPTHML(j)=ZOAV(j) end do ! j print *,' RCO2' - print 5001,((Rco2(j,k),j=1,jm),k=1,LMO) +! print 5001,((Rco2(j,k),j=1,jm),k=1,LMO) + print 5001,((Rco2(j,k)*1.e2,j=1,jm),k=1,LMO) dzog(1)=10./SQRT(1.7010587) dzg(2)=10. do l=2,lmo-1 @@ -335,8 +336,6 @@ READ (546) 245 continue endif - WRITE(503) OFFSSW 17.1 - REWIND 503 17.2 c CALL FRTR0(IO) 18. KBGN=KINC+1 18.5 KM2=KM*2-1 18.51 @@ -365,15 +364,16 @@ NSTEP2=NSTEP 29.6 MRCHT=0. 29.7 ITAU=(NSTEP+NSTEP0)*IDTHR 30. - TAU=FLOAT(ITAU)/XINT 31. +cjrs changed to dfloat 8/2/07 + TAU=DFLOAT(ITAU)/XINT 31. IDAY=1+ITAU/I24 32. TOFDAY=(ITAU-(IDAY-1)*I24)/XINT 33. - if(ISTART.eq.2.or.ISTRT1.eq.0.and..not.CONTRR)then - do 458 j=1,JM - TSURFD(j)=0. - TSURFT(j)=0. - 458 continue - endif +! if(ISTART.eq.2.or.ISTRT1.eq.0.and..not.CONTRR)then +! do 458 j=1,JM +! TSURFD(j)=0. +! TSURFT(j)=0. +! 458 continue +! endif if(JDATE.eq.100)then print *,JDATE,JMONTH,JYEAR print *,' main before daily0' @@ -388,8 +388,10 @@ endif CALL DAILY_NEW0 34. print *,' Main after DAILYNEW0 JYEAR=',JYEAR - print *,"DTSURF" - print *,DTSURF + print *,"DT2MGL" + print *,DT2MGL + print *,"DT2MLD" + print *,DT2MLD #if( !defined OCEAN_3D&& !defined ML_2D ) CALL DAILY_OCEAN print *,' AFTER DAILY_OCEAN IDAY=',IDAY,' IYEAR=',IYEAR @@ -438,25 +440,29 @@ JDATEATM=JDATE JYEARATM=JYEAR C +#if ( defined CPL_OCEANCO2 ) + do j=1,jm + fluxco2(j)=0.0 + enddo +#endif #if ( defined CPL_CHEM) && ( defined CPL_TEM ) C For TEM if(ISTRT1.eq.0) then c New run c Reading from flin_nep read(537)adupt,temco2 - else + else c Restart of the run c Reading from last_nep - read(367)adupt,temco2 -! & ,temch4,temn2o - rewind 367 - endif -! -! adupt= 1.459814341652516 -! adupt= 0.9078891180588442 -! adupt= 0.25 -! adupt= -0.1123070421398009 -! +cjrs file previously opened in input.F + read(876)adupt,temco2 +C CLOSE(876) + rewind 876 + endif + +cjrs next line per Andrei instruction 10/12/07 + adupt= 0.0 + aduptd=adupt/(365.*JM) temnepgl=0.0 do j=1,jm @@ -473,19 +479,20 @@ elseif(LMO.eq.12) then call ODIFS12 else - Print *,' Wromng LMO',LMO + Print *,' Wrong LMO',LMO stop endif endif #endif -#if (defined PREDICTED_GASES) +!#if (defined PREDICTED_GASES) #if (defined CPL_TEM || defined CPL_OCEANCO2 ) if(OBSFOR) then call obsco2(iyear,imontha) mnobco2=imonth endif #endif -#endif +!#endif +CJRS removed below from ocean_3d #ifdef ML_2D do j=1,jm do i=1,io @@ -524,8 +531,10 @@ enddo endif #endif +#if (defined CPL_TEM || defined CPL_OCEANCO2 ) print *,'ATM_CO2' print *,atm_co2 +#endif JDAYLAST=-1 ncallclm=0 NOCLM=.true. @@ -533,6 +542,7 @@ NOCLM=.false. #endif print *,' atmosphere DTATM=',DTATM + print *,' It is running' print *,'End of atmospheric model initialization' print *,' ' print *,' ' @@ -550,8 +560,6 @@ c HPRNT=TAU.ge.17520.0.and.TAU.lt.17545.0 c print *,' TAUE=',TAUE #if ( defined OCEAN_3D || defined ML_2D) -C print *,'TAU,DTATM,TAUE: ', TAU,DTATM,TAUE - CALL OCEAN4ATM do j=1,jm0 tauu(j)=0. tauv(j)=0. @@ -569,12 +577,13 @@ dhfidtgeq(j)=0. devidtgeq(j)=0. tempr(j)=0. +cjrs change var name to arunoff arunoff(j)=0. solarinc_ice(j)=0. solarnet_ice(j)=0. solarinc_ocean(j)=0. solarnet_ocean(j)=0. -Cjrs surfpr(j)=0. +Cjrs not used anymore (?) surfpr(j)=0. naveo(j)=0. navei(j)=0. navrad(j)=0. @@ -588,10 +597,23 @@ c enddo #endif +#ifdef OCEAN_3D +C get data from atm-ocean common block + do j=1,jm0 + ODATA(1,j,1)=mmsst(j) + ODATA(1,j,2)=mmfice(j) + GDATA(1,j,3)=mmtice(j) + GDATA(1,j,1)=mmsnowm(j) + ODATA(1,j,3)=mmicem(j) + GDATA(1,j,7)=0.5*(mmtice2(j)+mmtice1(j)) +# ifdef CPL_OCEANCO2 + fluxco2(j)=fluxco2(j) + dtatm*3600.*mmco2flux(j) +# endif + enddo +#endif WLMMAX=0.0 C 100 IF(.NOT.EVENT(TAUT)) GO TO 200 46. -C print *,' atmosphere TAU=',tau c HPRNT=TAU.ge.17520.00 NSTEP1=NSTEP 46.5 C**** WRITE RESTART INFORMATION ONTO DISK 47. @@ -614,7 +636,7 @@ if(TRANSR)then WRITE(KDISK0) AEXP,TAU,JC,C,RC,KEYNR,U,V,T,P,Q,ODATA,GDATA,BLDATA, 50. * RQT,SRHR,TRHR,(AJ(K,1),K=1,KACC),TAU,TSSFC,CKS,CKN,WMGE,TPRIM2 51. - * ,MRCHT,TRSURF,SRSURF,TSURFT,TSURFW,DWAV0 + * ,MRCHT,TRSURF,SRSURF,TLANDW,TSURFW,DWAV0 * ,TG3M,RTGO,STG3,DTG3 print *,' STG3' print 5001,(STG3(1,j),j=1,JM0) @@ -627,7 +649,7 @@ else WRITE(KDISK0) AEXP,TAU,JC,C,RC,KEYNR,U,V,T,P,Q,ODATA,GDATA,BLDATA, 50. * RQT,SRHR,TRHR,(AJ(K,1),K=1,KACC),TAU,TSSFC,CKS,CKN,WMGE,TPRIM2 51. - * ,MRCHT,TRSURF,SRSURF,TSURFT,TSURFW,DWAV0 + * ,MRCHT,TRSURF,SRSURF,TLANDW,TSURFW,DWAV0 endif REWIND KDISK0 52. end if ! ISTART.eq.2 @@ -637,13 +659,10 @@ PERCNT=100.*MELSE/(MSTART-MNOW+1.E-5) 56. MLAST=MNOW 59. C**** TEST FOR TERMINATION OF RUN 60. - 200 READ (503,END=210) LABSSW 61. + 200 continue c HPRNT=TAU.gt.45.0.and.TAU.lt.60.0 c HPRNT=TAU.gt.470.0.and.TAU.lt.550.0 NCOMP=0 - 210 REWIND 503 61.1 - IF(LABSSW.EQ.LABEL1) KSS6=1 61.2 - IF(KSS6.EQ.1) GO TO 800 62. IF(TAU+.06125.GE.TAUE) GO TO 820 63. JDAY00=JDAY C**** IF TIME TO ZERO OUT DIAGNOSTIC ACCUMULATING ARRAYS, DO SO 64. @@ -902,7 +921,8 @@ call chemmass66(1.0, 1.0,zco2,zco2mass) - call chemmass6(150.0,1.0,xn2o,xn2omass) +!call chemmass6(150.0,1.0,xn2o,xn2omass) + call chemmass6(120.0,1.0,xn2o,xn2omass) call chemmass2(1.0,ch4, ch4mass ) ! === if hfc, pfc, and sf6 are included: @@ -1030,7 +1050,6 @@ CALL RADIA endif #endif - if(HPRNT)then print *,' main after radia',' TAU=',TAU #include "PRNT.COM" @@ -1061,9 +1080,12 @@ print *,' main after surf4clm',' TAU=',TAU #include "PRNT.COM" endif + i=1 do j=1,jm - pcpl4clm(j)=pcpl4clm(j)*prlnd2total(j,mndriver) - pcpc4clm(j)=pcpc4clm(j)*prlnd2total(j,mndriver) + pcpl4clm(i,j)=pcpl4clm(i,j)*prlnd2total(j,mndriver) + & *3600./(NDYN*DT) + pcpc4clm(i,j)=pcpc4clm(i,j)*prlnd2total(j,mndriver) + & *3600./(NDYN*DT) enddo ! print *,' main after surf4clm',' TAU=',TAU ! print ('2(12f7.2,/,11f7.2,/)'),ps4clm,pcpl4clm, @@ -1100,12 +1122,12 @@ endif c if(JYEAR.gt.20)then -c write (934),tau,tsoiclm,snwdclm,snwcclm, -c & lwuclm,tref2mclm,tflxclm,tgndclm, -c & lhfclm,shfclm,tauxclm,tauyclm, -c & asdirclm,aldirclm,asdifclm,aldifclm, -c & sroclm,ssrclm,glrclm -c &,h2olclm,h2oiclm +! write (934),tau,tsoiclm,snwdclm,snwcclm, +! & lwuclm,tref2mclm,tflxclm,tgndclm, +! & lhfclm,shfclm,tauxclm,tauyclm, +! & asdirclm,aldirclm,asdifclm,aldifclm, +! & sroclm,ssrclm,glrclm +! &,h2olclm,h2oiclm c endif ! print *,' main after clm4mit2d',' TAU=',TAU ! print ('2(12f7.2,/,11f7.2,/)'),tsoiclm,snwdclm,snwcclm, @@ -1231,7 +1253,8 @@ C**** 189. 500 NSTEP=NSTEP+NDYN 190. ITAU=(NSTEP+NSTEP0)*IDTHR 191. - TAU=FLOAT(ITAU)/XINT 192. +cJRS fix to DFLOAT 8/2/07 + TAU=DFLOAT(ITAU)/XINT 192. IDAY=1+ITAU/I24 193. TOFDAYPR=TOFDAY+1.00 TOFDAY=(ITAU-(IDAY-1)*I24)/XINT 194. @@ -1239,6 +1262,7 @@ C 196. do J=1,JM0 TSURFW(J)=TSURFD(J) + TLANDW(J)=TLANDD(J) enddo JDATECLM=JDATE @@ -1363,6 +1387,9 @@ pC_atm(j)=zco2(1,j,1) & *28.97296245/44.0*1.e-9 !ppb(m) to kg per volume base + + atm_co2(j)=pC_atm(j)*1.e6 + enddo ! j c c ------- @@ -1401,32 +1428,26 @@ c print '12f7.1,/,2(11f7.1,/,),12f7.1',(pC_atm(j)*1.e6,j=1,jm) c print '12f7.1,/,2(11f7.1,/,),12f7.1',(rco2(j,1),j=1,jm) c ncallgary=ncallgary+1 - call carb_mxdlyr_chem(focean) - call carb_airsea_flx +! 10/28/06 +! call carb_mxdlyr_chem(focean) +! call carb_airsea_flx +! +! 3D ocean chemistry + call carb_chem_ocmip(focean) + call carb_airsea_flx(dtco2) +! 3D ocean chemistry +! 10/28/06 c print *,'FCO2 ncallgary=',ncallgary c print '12f7.1,/,2(11f7.1,/,),12f7.1', c & (fluxco2(j)*12.e-15*365.,j=1,jm) #endif -# if ( defined OCEAN_3D) - SECDAY=24.*3600. -c print *,'CO2F form ocean' -c print *,mmco2flux -Cjrs fluxco2(1)=SECDAY*mmco2flux(1) - fluxco2(1)=SECDAY*mmco2flux(2) - do j=2,jm-1 - fluxco2(j)=SECDAY*mmco2flux(j-1) - enddo - fluxco2(JM)=SECDAY*mmco2flux(JM-1) -Cjrs fluxco2(JM)=SECDAY*mmco2flux(JMOCEAN) -# endif - C For ocean carbon model c Annual oceanic CO2 uptake do j=1,jm OCUPT=OCUPT+fluxco2(j) enddo -c print *,' OCUPT=',OCUPT*12.e-15 +! print *,' OCUPT=',OCUPT*12.e-15 #if ( defined CPL_CHEM ) ! @@ -1440,6 +1461,9 @@ ! #endif + do j=1,jm + fluxco2(j)=0.0 + enddo #endif #if ( defined CPL_TEM ) @@ -1484,8 +1508,19 @@ #if ( defined CPL_OCEANCO2 && defined ML_2D ) C For OCM - dtco2=3600.*24. - call diffusco2(lmo,jm,dtco2,0.5,edzon,depthml,focean, +! dtco2=3600.*24. +! cfkvct=1.0 +! if (JYEAR.ge.1991) then begin +! if (JYEAR.le.2100) then begin +! cfkvct=(1.0*(2100-JYEAR)+0.25*(JYEAR-1990))/110. +! esle +! cfkvct=0.25 +! endif +! endif +! do j=1,jm +! edzcatr(j)=cfkvct*edzcar(j) +! enddo + call diffusco2(lmo,jm,dtco2,0.5,edzcart,depthml,focean, & dzg,dzog,rco2) call hdocean(rco2,focean,dxv,dyv,DXYP,depthml,edohd,dtco2) call avegary @@ -1724,8 +1759,8 @@ c do 5287 j=1,JM+3 c GBUDG(j,38,1)=GBUDG(j,37,1)*1013./GBUDG(jm+3,37,1) c5287 continue - print *,'FRMDICE' - print '6(1PE12.4)',FRMDICE +! print *,'FRMDICE' +! print '6(1PE12.4)',FRMDICE ENKE=0.0 ENPT=0.0 do ii=1,4 @@ -1776,13 +1811,13 @@ WRITE(KDISK0) AEXP,TAU,JC,C,RC,KEYNR,U,V,T,P,Q,ODATA,GDATA, * BLDATA, * RQT,SRHR,TRHR,(AJ(K,1),K=1,KACC),TAU,TSSFC,CKS,CKN,WMGE,TPRIM2 - * ,MRCHT,TRSURF,SRSURF,TSURFT,TSURFW,DWAV0 + * ,MRCHT,TRSURF,SRSURF,TLANDW,TSURFW,DWAV0 * ,TG3M,RTGO,STG3,DTG3 else WRITE(KDISK0) AEXP,TAU,JC,C,RC,KEYNR,U,V,T,P,Q,ODATA,GDATA, * BLDATA, * RQT,SRHR,TRHR,(AJ(K,1),K=1,KACC),TAU,TSSFC,CKS,CKN,WMGE,TPRIM2 - * ,MRCHT,TRSURF,SRSURF,TSURFT,TSURFW,DWAV0 + * ,MRCHT,TRSURF,SRSURF,TLANDW,TSURFW,DWAV0 c print *,' TSURFT' c print 5001,TSURFT c print *,' TSURFW' @@ -1803,7 +1838,7 @@ IF(TAU.GT.TAUX+3240.) GO TO 683 287. 685 WRITE (KCOPY) TAU,JC,C,RC,KEYNR,U,V,T,P,Q,ODATA,GDATA,BLDATA, 287.5 * RQT,SRHR,TRHR,(AJ(K,1),K=1,KACC),TAU,TSSFC,CKS,CKN 288. - * ,WMGE,TPRIM2,MRCHT,TRSURF,SRSURF,TSURFT,TSURFW,DWAV0 + * ,WMGE,TPRIM2,MRCHT,TRSURF,SRSURF,TLANDW,TSURFW,DWAV0 * ,TG3M,RTGO,STG3,DTG3 REWIND KCOPY 288.5 C 690 changed to 691 02/21/2003 @@ -1830,7 +1865,7 @@ REWIND KDISK0 303. WRITE(KDISK0) AEXP,TAU,JC,C,RC,KEYNR,U,V,T,P,Q,ODATA,GDATA,BLDATA, 304. * RQT,SRHR,TRHR,(AJ(K,1),K=1,KACC),TAU,TSSFC,CKS,CKN,WMGE,TPRIM2 305. - * ,MRCHT,TRSURF,SRSURF,TSURFT,TSURFW,DWAV0 + * ,MRCHT,TRSURF,SRSURF,TLANDW,TSURFW,DWAV0 end if C WRITE (6,908) 306. C**** RUN TERMINATED BECAUSE IT REACHED TAUE (OR SS6 WAS TURNED ON) 307. @@ -1841,12 +1876,19 @@ C DTATM time step of atm model in hours C precip and evap in mm/day or kg/m**2/day do j=1,jm0 -#if ( defined OCEAN_3D && defined CPL_OCEANCO2 ) - ncallatm=ncallatm+1 - co24ocean(j)=pC_atm(j)*1.e6 +Cjrs #if ( defined OCEAN_3D && defined CPL_OCEANCO2 ) +#ifdef OCEAN_3D +!jrs ncallatm=ncallatm+1 +! 020107 +! co24ocean(j)=pC_atm(j)*1.e6 +! jrs give CO2 even if ocn carbon off + co24ocean(j)=atm_co2(j) +# ifdef CPL_OCEANCO2 co24ocnan(j)=co24ocnan(j)+co24ocean(j) +# endif #endif #ifdef ML_2D +cjrs block only MD_2D rseaice(j)=ODATA(1,J,2) #endif tauu(j)=tauu(j)/(NSURF*DTATM) @@ -1858,7 +1900,7 @@ precip(j)=precip(j)*(1.-fland*prlnd2total(j,mndriver)) & /(1.-fland) endif -Cjrs surfpr(j)=surfpr(j)/(DTATM/24.) +Cjrs surfpr(j)=surfpr(j)/(DTATM/24.) c if(naveo(j).gt.0)then hfluxo(j)=NSURF*hfluxo(j)/(NDYN*DT*naveo(j)) @@ -1896,6 +1938,7 @@ endif c Runoff is a flux of water from land in mm/day c not for m**2 +cjrs change runoff to new name arunoff arunoff(j)=arunoff(j)/(DTATM/24.)*FDATA(1,j,2) & *DXYP(J) if(NWMGEA(J).gt.0)then @@ -1925,6 +1968,7 @@ CLAT=20.*TWOPI/360. do j=1,jm SLAND=SLAND+FDATA(1,j,2)*DXYP(J) +c jrs runoff->arunoff rungl=rungl+arunoff(j) if(LAT(J).lt.-CLAT)then runs=runs+arunoff(j) @@ -1938,6 +1982,7 @@ c print *,rungl/SLAND,rungl,runs,runt,runn c nmonth=JMNTH0 #ifdef ML_2D +c jrs only ML_2D nmonth=AMONTH(mndriver) #endif jdatefl=jdate-1 @@ -1986,6 +2031,7 @@ c & evai,hfluxo,dhfodtg,devodtg,hfluxi,dhfidtg,devidtg, c & solarinc_ice,solarnet_ice,rseaice #ifdef ML_2D +c jrs only ML_2D do j=1,jm osst(j)=ODATA(1,j,1) aoice(j)=ODATA(1,j,3) @@ -2008,6 +2054,7 @@ JYEARATM=JYEAR C #ifdef ML_2D +Cjrs change this block to only ML_2D IDAYM=IDAY JDAYM=JDAY JDATEM=JDATE @@ -2023,7 +2070,7 @@ c WRITE (6,905) TOFDAY,JDATE,JMONTH,JYEAR WRITE (6,905) TOFDAYPR,JDATEPR,JMONTHPR,JYEARPR endif - print *,'ncallclm=',ncallclm +cjrs print *,'ncallclm=',ncallclm JDAYLAST=JDAY c if(ncallclm.gt.6) stop c stop @@ -2031,7 +2078,6 @@ return C CALL ENQJOB 309. C CALL ENQJOB 310. - IF(KSS6.EQ.1) STOP 12 310.1 IF(IPFLAG.EQ.0) STOP 13 311. STOP 1 312. C**** 313.