! module temp use params,only: nlev use flds_hist,only: addfhist implicit none real,parameter :: ur=8.3144e+7 contains !----------------------------------------------------------------------- subroutine chemheat use flds_ratecoef use flds_atmos,only: xnn2,xno2,xno,xn4s,xnno,xn2d,xno3,xnoh,xnho2, | xnno2,rho,xnch4,xno1d,xnco2,xnh2,xnh2o,xno21d,xno21s,xnh2o2, | xnco,xnh,ch3o,ch2o,ch3ooh,ch3,ch3o2,hco,xnclo,xncl use flds_ionatm,only: xne,xiop2p,xiop2d,xinop,xinp,xiop,xin2p, | xio2p use cons,only: erg use flds_heat,only: qnnox,qno1d,qno21d,qno21s,qnhox,qnmth,qncl, | qnox,qnic,qn,qnnc implicit none ! ! Local: integer :: k real,parameter :: eps37 = 0.6 real :: xm,dum real :: qno1d_diag(nlev) ! diag for history file ! do k=1,nlev xm = xnn2(k)+xno2(k)+xno(k) qnnox(k) = (rb1(k)*xn4s(k)*xno2(k)*1.4+rb2(k)*1.84* | xn2d(k)*xno2(k)+rb3(k)*xn4s(k)*xnno(k)*3.25+ | rb4(k)*xn2d(k)*xno(k)*2.38+rb5(k)*xn2d(k)*xne(k)* | 2.38+rb6(k)*xn2d(k)*xnno(k)*5.63+rb8(k)*xn4s(k)* | xnoh(k)*2.10+rb9(k)*xnno(k)*xno3(k)*2.08+rb10(k)* | xnho2(k)*xnno(k)*0.35+rb11(k)*xno(k)*xnno2(k)*1.98+ | rb12(k)*xnno2(k)*xno3(k)*1.08+rb13(k)*xn4s(k)* | xnno2(k)*1.81)*erg/rho(k) ! ! qnnox(k) = 1.e-30 ! Roble experiment ! qno1d(k) = (rkm1(k)*xnn2(k)*1.96+rkm2a(k)*xno2(k)*0.33+ | rkm2b(k)*xno2(k)*1.96+rkm3(k)*xnh2o(k)*1.23+ | rkm4(k)*xnh2(k)*1.88+rkm5a(k)*xnch4(k)*1.85+ | rkm5b(k)*xnch4(k)*4.89+rkm6(k)*xnco2(k)*1.96+ | rkm7a(k)*xno3(k)*6.02+rkm7b(k)*xno3(k)*0.86)* | xno1d(k)*erg/rho(k) qno21d(k) = (rkm9(k)*xno2(k)*0.98+rkm10(k)*xnn2(k)*0.98+ | rkm11(k)*xno(k)*0.98)*xno21d(k)*erg/rho(k) qno21s(k) = (rkm13(k)*xnn2(k)*0.65+rkm14(k)*xnco2(k)*0.65+ | rkm15(k)*xno3(k)*0.65+rkm16(k)*xno(k)*0.65+ | rkm17(k)*xno2(k)*0.65)*xno21s(k)*erg/rho(k) qnhox(k) = (rkm25(k)*xno(k)*xnoh(k)*0.72+rkm26(k)*xno(k)* | xnho2(k)*2.33+rkm27(k)*xno(k)*xnh2o2(k)*3.44+ | rkm28(k)*xno(k)*xnh2(k)*0.08+ | rkm29(k)*xnoh(k)*xno3(k)*1.73+ | rkm30(k)*xnoh(k)*xnoh(k)*0.73+ | rkm31(k)*xnoh(k)*xnho2(k)*3.06+ | rkm32(k)*xnoh(k)*xnh2o2(k)*1.35+ | rkm33(k)*xnoh(k)*xnh2(k)*0.65+ | rkm34(k)*xnho2(k)*xno3(k)*1.23+ | rkm35(k)*xnho2(k)*xnho2(k)*1.71+ | rkm36(k)*xnh(k)*xno2(k)*xm*2.11+ | rkm37(k)*xnh(k)*xno3(k)*3.34*eps37+ | rkm38(k)*xnh(k)*xnho2(k)*2.41+ | rkm39(k)*xnh(k)*xnho2(k)*1.61+ | rkm40(k)*xnh(k)*xnho2(k)*2.34+ | rkm41(k)*xnh(k)*xnh(k)*xm*4.52+ | rkm42(k)*xnco(k)*xnoh(k)*1.07)*erg/rho(k) qnmth(k) = (gam1(k)*xnoh(k)*xnch4(k)*0.62+ | gam2a(k)*xno1d(k)*xnch4(k)*1.85+ | gam2b(k)*xno1d(k)*xnch4(k)*4.89+ | gam3(k)*ch3(k)*xno2(k)*xm*1.22+ | gam4(k)*xno(k)*ch3(k)*2.96+ | gam5(k)*xnno(k)*ch3o2(k)*0.72+ | gam6(k)*xnho2(k)*ch3o2(k)*1.74+ | gam7(k)*ch3o2(k)*ch3o2(k)*0.24+ | gam8(k)*xnoh(k)*ch3ooh(k)*1.32+ | gam9(k)*xno2(k)*ch3o(k)*1.15+ | gam10(k)*xnoh(k)*ch2o(k)*2.91+ | gam11(k)*xno(k)*ch2o(k)*2.18- | gam12(k)*xno2(k)*hco(k)*0.14- | gam13(k)*xno(k)*xnch4(k)*0.11+ | gam14(k)*xnco(k)*xno(k)*xm*5.51+ | gam15(k)*xno1d(k)*xnco2(k)*1.61)*erg/rho(k) qncl(k) = (del1(k)*xncl(k)*xno3(k)*1.68+ | del2(k)*xnclo(k)*xno(k)*2.38)*erg/rho(k) qnox(k) = (rkm20(k)*xno(k)*xno(k)*xm*5.12+ | rkm21(k)*xno(k)*xno2(k)*xno(k)*1.10+ | rkm22(k)*xno(k)*xno2(k)*xno2(k)*1.10+ | rkm23(k)*xno(k)*xno2(k)*xnn2(k)*1.10+ | rkm24(k)*xno(k)*xno3(k)*4.06)*erg/rho(k) qnic(k) = (rk1(k)*xno2(k)*xiop(k)*1.5555 | +rk2(k)*xnn2(k)*xiop(k)*1.0888 | +rk3(k)*xno(k)*xin2p(k)*0.7 | +rk4(k)*xn4s(k)*xio2p(k)*4.21 | +rk5(k)*xnno(k)*xio2p(k)*2.813 | +rk6(k)*xno2(k)*xinp(k)*2.486 | +rk7(k)*xno2(k)*xinp(k)*6.699 | +rk8(k)*xno(k)*xinp(k)*0.98 | +rk9(k)*xno2(k)*xin2p(k)*3.52 | +rk10(k)*xiop(k)*xn2d(k)*1.45 | +rk13(k)*xnco2(k)*xiop(k)*1.21 | +rk14(k)*xnh2(k)*xiop(k)*0.35 | +rk15(k)*xnh2o(k)*xiop(k)*1.04 | +rk16(k)*xnn2(k)*xiop2p(k)*3.02 | +rk17(k)*xnn2(k)*xiop2p(k)*0.70 | +rk18(k)*xno(k)*xiop2p(k)*5.0 | +rk19(k)*xiop2p(k)*xne(k)*5.0 | +rk20(k)*xiop2d(k)*xne(k)*1.69 | +rk23(k)*xno2(k)*xiop2p(k)*6.56 | +rk24(k)*xnn2(k)*xiop2d(k)*1.33 | +rk25(k)*xno(k)*xiop2d(k)*3.31 | +rk26(k)*xiop2d(k)*xne(k)*3.31 | +rk27(k)*xno2(k)*xiop2d(k)*4.865 | +alp1(k)*xinop(k)*xne(k)*0.854 | +alp2(k)*xio2p(k)*xne(k)*5.275 | +alp3(k)*xin2p(k)*xne(k)*3.678)*erg/rho(k) qnnc(k) = qnnox(k)+qno1d(k)+qno21d(k)+qno21s(k)+qnhox(k)+ | qnmth(k)+qncl(k)+qnox(k) qn(k)=qn(k)+qnnc(k)+qnic(k) enddo ! k=1,nlev ! call fprint8('chemheat',nlev, ! | qnnox , qno1d , qno21d , qno21s , qnhox , qnmth, qncl ,dum, ! | 'qnnox','qno1d','qno21d','qno21s','qnhox','qnmth','qncl' ,' ') ! call fprint8('chemheat',nlev, ! | qnox , qnic , qnnc , qn , dum , dum, dum ,dum, ! | 'qnox','qnic','qnnc','qn',' ',' ',' ' ,' ') ! ! Save diags to history output (sub addfhist in fhist.F): ! subroutine addfhist(f,varname,shortname,longname,units,nlevs, ! | itimedep,logplt) ! ! qno1d_diag = qno1d ! call addfhist(qno1d_diag,'QNO1D','QNO1D', ! | 'QNO1D DIAGNOSTIC','ergs/gm/s',nlev,1,1) ! call addfhist(qno1d_diag,'QNO1D_1d','QNO1D_1d', ! | 'QNO1D (1d)','ergs/gm/s',nlev,0,1) end subroutine chemheat !----------------------------------------------------------------------- subroutine glbmtemp(ut,dt) use flds_atmos,only: tn ! output use flds_modelz,only: zp use flds_heat,only: xco2c ! output use co2cooling,only: co2cool use flds_hist,only: addfhist implicit none ! ! Args: real,intent(in) :: ut,dt ! ! Local: real,parameter :: eps=1.e-2 integer,parameter :: npde=1 integer,parameter :: nwkcn=(24+6)*(nlev+2)+5*(nlev+2)+4, niwcn=4 integer :: iwcn(niwcn) real :: wkcn(nwkcn),zpts(nlev),usoln(npde,nlev),ddt,tm,tend real :: xlco2e(nlev),xlco2d(nlev) integer :: k,index,mblft(npde),mbrht(npde) integer,parameter :: itype=2, nbug=1 real :: xco2c_diag(nlev) ! diag only ! ! Common /syspar/ is retained from old code for pdedif solver: real :: uround integer :: min,mout common/syspar/ uround,min,mout ! mblft(1) = 1 mbrht(1) = 2 tm = 0. uround = 7.e-15 ! /syspar/ for pdedif min = 5 ! /syspar/ for pdedif mout = 6 ! /syspar/ for pdedif ddt = 1.e-3 index = 1 tend = dt ! ! Normal is tn(1)=225. ! Present day temperatue at 30 mb ! tn(1)=225. ! co2 doubled yr=2000 temperature at 30 mb tn(1)=215. ! old model (also used 225., 235.) ! co2 halved ice age temperature at 30 mb ! tn(1)=235. ! zpts(:) = zp(:) usoln(1,:) = tn(:) ! ! Fomichev co2 cooling is in co2cool.F. call co2cool(xlco2e,xlco2d) ! ! Set xco2c (flds_heat): do k=1,nlev xco2c(k) = -xlco2e(k) enddo ! k=1,nlev ! ! Save xco2c on output history file: ! subroutine addfhist(f,varname,shortname,longname,units,nlevs, ! | itimedep,logplt) ! xco2c_diag = xco2c xco2c_diag(1) = xco2c_diag(2) ! change lbc (which looks wrong in xco2c) ! call addfhist(xco2c_diag,'CO2COOL','CO2COOL','CO2 COOLING',' ', ! | nlev,1,1) ! ! PDE solver: call pdedif(func,bnbdy,npde,nlev,tend,eps,zpts,itype,mblft, | mbrht,nbug,tm,ddt,index,usoln,wkcn,iwcn,nwkcn,niwcn) tn(:) = usoln(1,:) ! ! tn diff with old code at tn(66) by 0.1 degrees. ! write(6,"('glbmtemp after pdedif: tm=',e12.4,' ddt=',e12.4, ! | ' tn=',/,(6e12.4))") tm,ddt,tn ! write(6,"('glbmtemp after pdedif: ut=',f5.2,' tn=',/,(6e12.4))") ! | ut,tn end subroutine glbmtemp !----------------------------------------------------------------------- subroutine func(npde,nzp2,zn,tm,un,uzn,an,bn,cn) ! ! User-provided subroutine for pdedif solver (see pdedif.F). ! use flds_atmos,only: tn,xnn2,xno2,xno,xnh,xnhe,rho,sht use flds_heat,only: | cmu,cp,ck,hcm,hce,gwheat,xircool,xmetac,xcool, ! output | qn,xco2c,xnoc,xo3pc,xno21dc,xno21sc ! input use flds_modelz,only: gz,zp use flds_fields,only: edy implicit none ! ! Args: integer,intent(in) :: npde,nzp2 real,intent(in) :: zn(nzp2),un(npde,nzp2),uzn(npde,nzp2),tm real,intent(out) :: an(npde,nzp2),bn(npde,nzp2),cn(npde,nzp2) ! ! Local: integer :: k,kk integer :: nlevp1=nlev+1, nlevp2=nlev+2, nlevm1=nlev-1 real,dimension(nlev) :: ac,edt,dedy,dsht,dck,dcp,drho, | dgz,fmt,fet real :: prso,dzit,dziu,xnt,tt,pn2,po,po2,gweff,gwcof,gwht1, | gwht2,dum ! ! Alternate values (commented out in old code) for prandl number ! prntl include 1., 3., 10., 1.e+20 ! real,parameter :: prntl = 1.0 ! prso=5.E-4 dzit=0.5 dziu=0.25 do kk=2,nlevp1 ! was do 1 k=kk-1 ! ! Set tn (module flds_atmos): tn(k)=un(1,kk) tt=un(1,kk)**0.69 xnt=xnn2(k)+xno2(k)+xno(k) ! ! Set rho (module flds_atmos): rho(k)=(28.*xnn2(k)+32.*xno2(k)+16.*xno(k)+xnh(k)+xnhe(k))* | 1.66E-24 pn2=xnn2(k)/xnt po=xno(k)/xnt po2=1.-pn2-po ! ! Set cp and ck (module flds_heat) ! cp units: ergs/(k-cm**3) ! ck units: ergs/(cm-k-sec) ! ck(k)=((po2+pn2)*56.+po*75.9)*tt cp(k)=0.5*ur*(po2*7./32.+pn2*7./28.+po*5./16.) ! cmu(k)=(po2*4.03+pn2*3.42+po*3.9)*1.e-6*tt ac(k)=gz(k)*exp(zp(k))/(prso*cp(k)) edt(k)=edy(k)/prntl ! ! an(npde,nzp2) output to pde solver: an(1,kk)=ac(k)*(ck(k)/sht(k)+edt(k)*sht(k)*cp(k)*rho(k)) enddo ! kk=2,nlevp1 ! do k=2,nlevm1 ! was do 2 dedy(k) = (edt(k+1)-edt(k-1))/dzit dsht(k) = (sht(k+1)-sht(k-1))/dzit dck(k) = (ck (k+1)-ck (k-1))/dzit dcp(k) = (cp (k+1)-cp (k-1))/dzit drho(k) = (rho(k+1)-rho(k-1))/dzit dgz(k) = (gz (k+1)-gz (k-1))/dzit enddo ! k=2,nlevm1 ! ! Lower boundary: dedy(1) = (edt(2)-edt(1))/dziu dck (1) = (ck (2)-ck (1))/dziu dsht(1) = (sht(2)-sht(1))/dziu dcp (1) = (cp (2)-cp (1))/dziu drho(1) = (rho(2)-rho(1))/dziu dgz (1) = (gz (2)-gz( 1))/dziu ! ! Upper boundary: dedy(nlev) = (edt(nlev)-edt(nlevm1))/dziu dck (nlev) = (ck (nlev)-ck (nlevm1))/dziu dsht(nlev) = (sht(nlev)-sht(nlevm1))/dziu dcp (nlev) = (cp (nlev)-cp (nlevm1))/dziu drho(nlev) = (rho(nlev)-rho(nlevm1))/dziu dgz (nlev) = (gz (nlev)-gz (nlevm1))/dziu ! do k=1,nlev ! was do 4 kk=k+1 fmt(k)=ck(k)*uzn(1,kk)/sht(k) fet(k)=edt(k)*sht(k)*sht(k)*cp(k)*rho(k)* | (gz(k)/cp(k)+uzn(1,kk)/sht(k)) enddo ! k=1,nlev ! ! Set hcm, hce (module flds_heat): do k=2,nlevm1 ! was do 5 hcm(k)=ac(k)*cp(k)*(fmt(k+1)-fmt(k-1))/dzit hce(k)=ac(k)*cp(k)*(fet(k+1)-fet(k-1))/dzit enddo ! ! hcm, hce lower boundary: hcm(1)=ac(1)*cp(1)*(fmt(2)-fmt(1))/dziu hce(1)=ac(1)*cp(1)*(fet(2)-fet(1))/dziu ! ! hcm,hce upper boundary: hcm(nlev)=ac(nlev)*cp(nlev)*(fmt(nlev)-fmt(nlevm1))/dziu hce(nlev)=ac(nlev)*cp(nlev)*(fet(nlev)-fet(nlevm1))/dziu ! gweff = 0.1 ! not used? (was in old code) ! ! Infra-red cooling: call ircool ! ! gwheat: gravity-wave heating (flds_heat): do kk=2,nlevp1 k=kk-1 gwcof=1.e-20 gwht1=gwcof*gz(k)/(sht(k)*un(1,kk)) gwht2=gwcof*gz(k)*gz(k)/(cp(k)*un(1,kk)) gwheat(k)=(gwht1+gwht2)*cp(k) ! ! bn(npde,nzp2) output to pde solver: bn(1,kk)=ac(k)*(dck(k)/sht(k)-ck(k)*dsht(k)/(sht(k)*sht(k))+ | sht(k)*cp(k)*rho(k)*dedy(k)+edt(k)*cp(k)*rho(k)*dsht(k)+ | edt(k)*sht(k)*rho(k)*dcp(k)+edt(k)*sht(k)*cp(k)*drho(k))+gwht1 xircool(k)=xco2c(k)+xnoc(k)+xo3pc(k) xmetac(k)=xno21dc(k)+xno21sc(k) xcool(k)=xircool(k)+xmetac(k) ! ! cn(npde,nzp2) output to pde solver: cn(1,kk)=(qn(k)-xcool(k))/cp(k)+(sht(k)*sht(k)*rho(k)*gz(k)* | dedy(k)+edt(k)*rho(k)*gz(k)*2.*sht(k)*dsht(k)+edt(k)*sht(k)* | sht(k)*gz(k)*drho(k)+edt(k)*sht(k)*sht(k)*rho(k)*dgz(k))*ac(k) | +gwht2 enddo ! kk=2,nlevp1 ! ! Boundaries of an,bn,cn output: an(1,1)=an(1,2) bn(1,1)=bn(1,2) cn(1,1)=cn(1,2) an(1,nzp2)=an(1,nzp2-1) bn(1,nzp2)=bn(1,nzp2-1) cn(1,nzp2)=cn(1,nzp2-1) ! ! call fprint8('func for pdedif',nlev, ! | an , bn , cn , dum , dum , dum, dum ,dum, ! | 'an','bn','cn',' ',' ',' ',' ' ,' ') end subroutine func !----------------------------------------------------------------------- subroutine ircool ! ! Infra-red cooling. Define xnoc, xo3pc (flds_heat) ! use flds_atmos,only: tn,xno,rho,xnno use flds_heat,only: xnoc,xo3pc ! output use cons,only: boltz implicit none integer :: k real :: xfac(nlev) real,parameter :: a10=13.3, hvno=3.726e-13, rk10=4.2e-11 ! xfac(1:29) = .01 xfac(30:nlev-53) = | (/0.025, 0.05, 0.075, 0.1 , 0.15, 0.2 , 0.3 , | 0.4 , 0.48, 0.55 , 0.64, 0.7 , 0.725, 0.75, 0.775/) xfac(nlev-53+1:nlev) = 0.8 ! do k=1,nlev xnoc(k)=hvno*xnno(k)*a10*((rk10*xno(k))/(rk10*xno(k)+a10))* | exp(-hvno/(boltz*tn(k)))/rho(k) ! xnoc(k) = 1.e-30 ! Roble experiment ! xo3pc(k)=(1.67e-18*xno(k)*exp(-228./tn(k))/(1.+0.6*exp(-228./ | tn(k))+0.2*exp(-325./tn(k))))/rho(k)*xfac(k) enddo ! k=1,nlev ! end subroutine ircool !----------------------------------------------------------------------- subroutine bnbdy(npde,tm,zl,zr,al,ar,bl,br) ! ! Called by pde solver in pdedif.F ! use flds_atmos,only: tn implicit none integer,intent(in) :: npde real,intent(in) :: tm,zl,zr real,intent(in) :: al(npde),ar(npde) real,intent(out) :: bl(npde),br(npde) ! bl(1)=tn(1) br(1)=0. end subroutine bnbdy !----------------------------------------------------------------------- subroutine rhozmod ! ! Reset rho, sht, and zpht. Called after glbmtemp. ! use flds_atmos,only: tn,xno,xno2,xnn2,xno,xnh,xnhe,rho,sht use flds_modelz,only: gz,zpht,zp use cons,only: rmn2,rmo2,rmo,rmh,rmhe implicit none integer :: k real :: xnt,dzp,xmb real,parameter :: eee=1.66026E-24 real,parameter :: boltz=1.38066e-16 real,parameter :: | wtmol_n2 = rmn2*eee, ! was wtmol(1) | wtmol_o2 = rmo2*eee, ! was wtmol(2) | wtmol_o = rmo *eee, ! was wtmol(3) | wtmol_he = rmhe*eee, | wtmol_h = rmh *eee ! do k=1,nlev ! was do 1 xnt=xnn2(k)+xno2(k)+xno(k) ! ! Reset rho in module flds_atmos: rho(k) = wtmol_n2*xnn2(k)+wtmol_o2*xno2(k)+wtmol_o*xno(k)+ | wtmol_h *xnh(k) +wtmol_he*xnhe(k) xmb=rho(k)/xnt ! ! Reset sht in module flds_atmos: sht(k)=boltz*tn(k)/(gz(k)*xmb) ! write(6,"('rhozmod: k=',i3,' tn=',e12.4,' gz=',e12.4,' xmb=', ! | e12.4,' boltz=',e12.4,' sht=',e12.4)") k,tn(k),gz(k),xmb, ! | boltz,sht(k) enddo ! k=1,nlev dzp = zp(2)-zp(1) ! dzp is local (was in modelp.h in old code) ! write(6,"('rhozmod: sht=',/,(6e12.4))") sht do k=2,nlev zpht(k)=zpht(k-1)+(sht(k)+sht(k-1))*0.5*dzp*1.e-5 enddo ! k=2,nlev ! write(6,"('rhozmod: rho=',/,(6e12.4))") rho ! write(6,"('rhozmod: zpht=',/,(6e12.4))") zpht end subroutine rhozmod !----------------------------------------------------------------------- end module temp