c ------------------------------------------------------------ subroutine tr_integrate(npt,nz,nstep,nstart,nzi,aarea, & trmf,hm,trint,trfrst,vollev,trlev,sqrlev,itrac_cons) c ------------------------------------------------------------ c subroutine to integrate tracer amount over model domain c taking into account that area varies at corners and edges; c c trmf: tracer concentration [same as tr()] c trint: integrated tracer content c [analog of en(npten,1)] c c Tracer integrated over 2-D planes: {trint <- tramt} c c trint(m,1,k)= tracer amount for k vertical layer c for tracer m c k: 1<->nz c trint(m,2,k)= tracer variance for k vertical level c for tracer m c sum(h*area*tr^2) c trint(m,3,k)= volume of the k'th level c c Tracer integerated over 3-D model space: c c trint(m,1,nz+1) = tracer amount over 3-d domain space c trint(m,2,nz+1) = tracer variance c trint(m,3,nz+1) = normalized volume for entire domain; c implicit real(a-h,o-z),integer(i-n) c include 'comm_tracer.h' include 'comm_para.h' include 'comm_new.h' dimension trmf(npt,nz),hm(npt,nz),aarea(npt,nz) dimension trint(ntrac,3,nz+1),trfrst(ntrac,3,nz+1) dimension vollev(nz), trlev(nz), sqrlev(nz) dimension nzi(npt) logical FIRST save FIRST data FIRST /.true./ do m=1,ntrac trtot = 0.0 sqrtot = 0.0 do k=1,nz trlev(k)=0.0 sqrlev(k)=0.0 vollev(k)=0.0 enddo do i=1,npt do k=1,nzi(i) volbox = aarea(i,k)*hm(i,k) trbox = trmf(i,k)*volbox sqrbox = trmf(i,k)*trbox vollev(k)= vollev(k)+ volbox trlev(k) = trlev(k) + trbox sqrlev(k)= sqrlev(k)+sqrbox sqrtot = sqrtot + sqrbox voltot = voltot + volbox trtot = trtot + trbox enddo enddo do k=1,nz trint(m,1,k)=trlev(k)/vollev(k) trint(m,2,k)=(sqrlev(k)/vollev(k))**0.5 enddo trint(m,1,nz+1)=trtot/voltot trint(m,2,nz+1)=(sqrtot/voltot)**0.5 enddo c write(*,*) trint(1,1,nz+1),trint(1,2,nz+1) cc Note: Be careful when trfrst=0!!! cc What happens for the bottom level when cc using sigma coordinates? cc if (FIRST) then do k=1,nz+1 do m=1,ntrac trfrst(m,1,k) = trint(m,1,k) trfrst(m,2,k) = trint(m,2,k) trfrst(m,3,k) = trint(m,3,k) enddo enddo else do m=1,ntrac do k=1,nz trint(m,1,k) = trint(m,1,k)/trfrst(m,1,k) trint(m,2,k) = trint(m,2,k)/trfrst(m,2,k) trint(m,3,k) = trint(m,3,k)/trfrst(m,3,k) enddo trint(m,1,nz+1) = trint(m,1,nz+1)/trfrst(m,1,nz+1) trint(m,2,nz+1) = trint(m,2,nz+1)/trfrst(m,2,nz+1) trint(m,3,nz+1) = trint(m,3,nz+1)/trfrst(m,3,nz+1) enddo endif if (.not. FIRST .and. itrac_cons.eq.1) then write(92,*) trint(1,1,nz+1) - 1. endif c if (FIRST) then c do i=1,500 c write(*,*) aarea(i,1), i c enddo c endif if (FIRST) then FIRST = .not. FIRST endif return end c------------------------------------------------------- subroutine rstrt_test(npt,nz,ntrac,irest,tr,trm) c------------------------------------------------------- real tr(npt,nz,ntrac),trm(npt,nz,ntrac) write(*,*) " " write(*,*) " w/i rstrt_test: " write(*,*) " irest= ", irest write(*,*) " tr(npt,1,1)= ", tr(npt,1,1) c write(*,*) " tr(npt,1,2)= ", tr(npt,1,2) write(*,*) " " return end c-------------------------------------------------- subroutine hflx_pert(npt,nz,nx,ny,nstep,yy) c-------------------------------------------------- implicit real(a-h,o-z),integer(i-n) include 'comm_data.h' include 'comm_new.h' include 'comm_tracer.h' dimension yy(ny) qcon_inv = 1./QCON add_term = qcon_inv*hfprt c apply perturbation (hfprt) c uniformly over domain; if (ihfprt .eq. 1) then do i=1,npt q(i) = q(i) + add_term enddo c apply perturbation poleward c of specified latitude (hfprt_lat) c but w/ hemispheric symmetry elseif (ihfprt .eq. 2) then do i=1,npt j=(iox(i)-1)/nx + 1 if (abs(yy(j)).ge.hfprt_lat) then q(i) = q(i) + add_term endif enddo c apply only in northern hemisphere elseif (ihfprt .eq. 3) then do i=1,npt j=(iox(i)-1)/nx + 1 if (yy(j).ge.hfprt_lat) then q(i) = q(i) + add_term endif enddo c apply only in southern hemisphere elseif (ihfprt .eq. 4) then do i=1,npt j=(iox(i)-1)/nx + 1 if (yy(j).le.hfprt_lat) then q(i) = q(i) + add_term endif enddo endif return end c-------------------------------------------------- subroutine force_tritium2(npt,nz,ntrac,nstep,nxp,nyp,dt, & juljar,rjuljar,nzi, & yy,t,evp,prc,rlh,abw,flx,tr,ftr,hm,iox,tpw,tpf) c-------------------------------------------------- implicit real(a-h,o-z),integer(i-n) include 'comm_tracer.h' parameter(spyr=365.*86400.) parameter(ss1972=1958.6/293.8) dimension yy(1) dimension evp(npt,1),prc(npt,1),rlh(npt,1),abw(npt,1) dimension flx(1) dimension t(npt,nz),tr(npt,nz,ntrac),ftr(npt,nz,ntrac),hm(npt,nz) dimension iox(npt),nzi(npt) dimension tpw(npt,3),tpf(npt,1) pi = asin(1.)*2. rlambda = 12.43*365.*24.*3600. decay_term = exp(-dt*alog(2.)/rlambda) - 1. decay_term = decay_term/dt c RADIOACTIVE DECAY do i=1,npt do k=1,nzi(i) trit_sink = -hm(i,k)*tr(i,k,1)*decay_term ftr(i,k,1) = ftr(i,k,1)-trit_sink ftr(i,k,2) = ftr(i,k,2)+trit_sink enddo enddo c AIR-SEA EXCHANGE c from Christoph's subroutine c "gasexc2.F" c Follows equations in Wanninkhof c JGR 1992 do i=1,npt tcel = t(i,1) tcel2 = tcel*tcel tcel3 = tcel*tcel2 asea = 410.14 bsea = 20.503 csea = 0.53175 dsea = 0.0060111 scsea = asea-bsea*tcel+csea*tcel2-dsea*tcel3 c...eq (1) in wanninkhof (1992) cm/hr rk(i) = 0.39*tpw(i,1)**2/sqrt(scsea/660.) c rk(i) = 0.39*5.**2/sqrt(scsea/660.) c rk(i) = 0.39*abswin(i)**2/sqrt(scsea/660.) c...convert to m/second rk(i) = rk(i)/(100.*3600.) term = rk(i)*tr(i,1,2) ftr(i,1,2) = ftr(i,1,2) -term enddo return end c-------------------------------------------------- subroutine force_tritium(npt,nz,ntrac,nstep,nxp,nyp,dt, & juljar,rjuljar,nzi, & yy,t,evp,prc,rlh,abw,flx,tr,ftr,hm,iox,tpw,tpf) c-------------------------------------------------- implicit real(a-h,o-z),integer(i-n) include 'comm_tracer.h' parameter(spyr=365.*86400.) parameter(ss1972=1958.6/293.8) parameter(trit_offset=100.0) dimension yy(1) dimension evp(npt,1),prc(npt,1),rlh(npt,1),abw(npt,1) dimension flx(1) dimension t(npt,1),tr(npt,1,1),ftr(npt,1,1),hm(npt,1) dimension iox(npt),nzi(npt) dimension tpw(npt,1),tpf(npt,1) pi = asin(1.)*2. rlambda = 12.43*365.*24.*3600. decay_term = exp(-dt*alog(2.)/rlambda) - 1. decay_term = decay_term/dt c Find year in W&R data which corresponds c to model year c [note: eventually I should interpolate c between W&R annual mean values] if ((iforc_tr.eq.61 .or. iforc_tr.eq.62).or. & (iforc_tr.eq.63 .and. juljar.lt.1960)) then nin=0 do n=1,newr if (juljar .eq. nint(souryr(n))) nin=n enddo endif c i/o stuff for evap,precip,relhum if (iforc_tr.eq.62 .or. iforc_tr.eq.63) then call it_catch(nevap,tdoney,nstep,it1,it2,tscl) if (it2.ne.idoney) then idoney=it2 do i=1,npt prc(i,1)=prc(i,2) evp(i,1)=evp(i,2) rlh(i,1)=rlh(i,2) abw(i,1)=abw(i,2) enddo call read_zt(idf_donprcp,lprecip,npt,1,it2,'precip', & prc(i,2)) call read_zt(idf_donevp,levap,npt,1,it2,'evap', & evp(i,2)) call read_zt(idf_donrh,lrhum,npt,1,it2,'relhum', & rlh(i,2)) call read_zt(idf_donabwn,labwn,npt,1,it2,'abswin', & abw(i,2)) endif do i=1,npt prc(i,3) = prc(i,1) + tscl*(prc(i,2)-prc(i,1)) evp(i,3) = evp(i,1) + tscl*(evp(i,2)-evp(i,1)) rlh(i,3) = rlh(i,1) + tscl*(rlh(i,2)-rlh(i,1)) abw(i,3) = abw(i,1) + tscl*(abw(i,2)-abw(i,1)) enddo endif if (iforc_tr.eq.61) then c Find latitude index "jin" from W&R data c which corresponds to model's "i" c gridpoint; c [note: eventally do spatial interp.] do i=1,npt j = ((iox(i)-1)/nxp)+1 do jj=2,jewr-1 phitop = souphi(jj) + 2.5 phibot = souphi(jj) - 2.5 if (yy(j).le.phitop .and. yy(j).gt.phibot) jin=jj enddo avsis = 0.0 if (sousp(jin,2) .gt. 0.0) then avsis50 = sousp(jin,2) endif if (souphi(jin) .ge. 0.0) then cp(i) = cp50n(nin)*avsis50 else cp(i) = cp50s(nin)*avsis50*ss1972 endif hrel=0.74 fac1 = (hrel/(1.-hrel))/fralpha fac2 = 1./(fralpha*(1.-hrel)) depni = (soup(jin,2) + fac1*soue(jin,2))*cp(i) dvni = (sourv(jin,2)/soua(jin,2))*3.*cp(i) third_term = soue(jin,2)*fac2*(tr(i,1,1)) c third_term = soue(jin,2)*fac2*(trit_offset-tr(i,1,1)) source_term = depni + dvni - third_term flx(i) = source_term/hm(i,1) c flx(i) = 36.0/hm(i,1) ftr(i,1,1) = ftr(i,1,1) + flx(i)*hm(i,1)/(86400.*365.) c ftr(i,1,1) = ftr(i,1,1) + source_term/(86400.*365.) enddo endif if (iforc_tr.eq.62 .or. (iforc_tr.eq.63 .and. & juljar.lt.1960)) then do i=1,npt j = ((iox(i)-1)/nxp)+1 do jj=2,jewr-1 phitop = souphi(jj) + 2.5 phibot = souphi(jj) - 2.5 if (yy(j).le.phitop .and. yy(j).gt.phibot) jin=jj enddo avsis = 0.0 if (sousp(jin,2) .gt. 0.0) then avsis50 = sousp(jin,2) endif if (souphi(jin) .ge. 0.0) then cp(i) = cp50n(nin)*avsis50 else cp(i) = cp50s(nin)*avsis50*ss1972 endif hrel=0.74 term1 = prc(i,3)*cp(i) term2 = evp(i,3)*rlh(i,3)*cp(i)/(fralpha*(1.-rlh(i,3))) term3 = evp(i,3)*(tr(i,1,1))/(fralpha*(1.-rlh(i,3))) c term3 = evp(i,3)*(tr(i,1,1)-trit_offset)/(fralpha*(1.-rlh(i,3))) source_term = term1 + term2 - term3 trtflx1(i) = term1/hm(i,1) trtflx2(i) = term2/hm(i,1) trtflx3(i) = term3/hm(i,1) flx(i) = trtflx1(i) + trtflx2(i) - trtflx3(i) c flx(i) = 36.0/hm(i,1) ftr(i,1,1) = ftr(i,1,1) + hm(i,1)*flx(i)/(86400.*365.) c ftr(i,1,1) = ftr(i,1,1) + source_term/(86400.*365.) enddo endif if (iforc_tr.eq.63 .and. juljar.ge.1960) then ndonin = 0 do n=1,nedon idonyr = nint(donyr(n)-0.5) if (juljar .eq. idonyr) ndonin = n enddo decyr = rjuljar - real(juljar) do i=1,npt cp(i)=donf1(i)*cptdon(ndonin,1)+donf2(i)*cptdon(ndonin,2) cp(i)=cp(i)*(1.+donam(i)*cos(2.*pi*(decyr-donphz(i)))) term1 = prc(i,3)*cp(i) term2 = evp(i,3)*cp(i)*rlh(i,3)/(fralpha*(1.-rlh(i,3))) term3 = evp(i,3)*(tr(i,1,1))/(fralpha*(1.-rlh(i,3))) c term3 = evp(i,3)*(tr(i,1,1)-trit_offset)/(fralpha*(1.-rlh(i,3))) source_term = term1 + term2 - term3 trtflx1(i) = term1/hm(i,1) trtflx2(i) = term2/hm(i,1) trtflx3(i) = term3/hm(i,1) flx(i) = trtflx1(i) + trtflx2(i) - trtflx3(i) c flx(i) = 36.0/hm(i,1) ftr(i,1,1) = ftr(i,1,1) + hm(i,1)*flx(i)/(86400.*365.) c ftr(i,1,1) = ftr(i,1,1) + source_term/(86400.*365.) enddo endif c AIR-SEA EXCHANGE c from Christoph's subroutine c "gasexc2.F" c Follows equations in Wanninkhof c JGR 1992 if (igas_ex .ge. 3) then do i=1,npt tcel = t(i,1) tcel2 = tcel*tcel tcel3 = tcel*tcel2 asea = 410.14 bsea = 20.503 csea = 0.53175 dsea = 0.0060111 scsea = asea-bsea*tcel+csea*tcel2-dsea*tcel3 c...eq (1) in wanninkhof (1992) cm/hr rk(i) = 0.39*tpw(i,1)**2/sqrt(scsea/660.) c if (i.eq.1000) then c write(*,*) " " c write(*,*) " tcel= ", tcel c write(*,*) " wnsp= ", tpw(i,1) c write(*,*) " rk= ", rk(i) c write(*,*) " rk/(100.*3600.)= ", rk(i)/(100.*3600.) c write(*,*) " [piston velocity in m/s]" c write(*,*) " rk*24/100.= ", rk(i)*24./100. c write(*,*) " [piston velocity in m/day]" c write(*,*) " scsea= ", scsea c endif c rk(i) = 0.39*5.**2/sqrt(scsea/660.) c rk(i) = 0.39*abswin(i)**2/sqrt(scsea/660.) c...convert to m/second rk(i) = rk(i)/(100.*3600.) term = rk(i)*tr(i,1,2) ftr(i,1,2) = ftr(i,1,2) -term enddo endif c RADIOACTIVE DECAY do i=1,npt do k=1,nzi(i) trit_sink = -hm(i,k)*tr(i,k,1)*decay_term ftr(i,k,1) = ftr(i,k,1)-trit_sink ftr(i,k,2) = ftr(i,k,2)+trit_sink c dcdt = decay_term*tr(i,k,1)/dt c term2 = hm(i,k)*dcdt c ftr(i,k,1) = ftr(i,k,1)+term2 c ftr(i,k,2) = ftr(i,k,2)-term2 enddo enddo return end c-------------------------------------------------- subroutine force2_tritium(npt,nz,ntrac,nstep,nzi,ftr,tr,hm,flx) c-------------------------------------------------- implicit real(a-h,o-z),integer(i-n) include 'comm_tracer.h' dimension nzi(npt) dimension flx(1) dimension tr(npt,nz,ntrac),ftr(npt,nz,ntrac),hm(npt,nz) c for the following example, c fake=3.0 c (units: TU m/yr) c means that a ten meter layer will c increase its tritium concentration c by 0.3 units over one year c fake_flux = 36.5 do i=1,npt flx(i) = fake_flux ftr(i,1,1) = ftr(i,1,1) + flx(i)/(86400.*365.) enddo c do i=1,npt c flx(i) = 36.0 c do k=1,nzi(i) c ftr(i,k,1) = ftr(i,k,1) + hm(i,1)*flx(i)/(86400.*365.) c enddo c enddo c do i=1,npt c flx(i) = 36.0 c ftr(i,1,1) = ftr(i,1,1) + hm(i,1)*flx(i)/(86400.*365.) c enddo return end c------------------------------------------------------------ subroutine force_tracer(tenso,npt,nz,ntrac,nstep,nx,ny,iv_bot, & rjuljar,juljar,dnt, & nzi,tr,ftr,t,h,tpw, & sal,pdens,yy,iox,tpf) c------------------------------------------------------------ implicit real(a-h,o-z),integer(i-n) include 'comm_tracer.h' include 'comm_pbl.h' dimension tpw(npt,3) dimension tr(npt,nz,ntrac),ftr(npt,nz,ntrac) dimension h(npt,nz),t(npt,nz),sal(npt,nz),yy(ny) dimension pdens(npt,nz) dimension nzi(npt) dimension iox(npt) dimension flux(100),depth(100),term(100) dimension tpf(1) parameter (factor = 0.5) parameter (secpmonth = 30.*24.*3600.) parameter (tau = secpmonth) parameter (tau_year = 1./(24.*365.*3600.)) c parameter (factor_c14 = 1./(24.*3600.*365.*7.5)) parameter (factor_age = 1./(24.*3600.*365.)) parameter (secpyr = 24.*3600.*365.) parameter (trstar = 1.0) parameter (factor_10day = 1/(10.*86400.)) parameter (factor_1month = 1/((365./12.)*86400.)) c Converting to indices for X/Y position: c j = ((iox(k)-1)/nx) + 1 c i = iox(k) - (j-1)*nx rjuljar = tenso/12.+1960. juljar = int(rjuljar) c for other options, see tracer.f and trac_init.f do itrac = 1, ntrac itracp = itrac + 1 iforc_tr = forctr(itrac) if (iforc_tr.eq.0) goto 100 c WINDSPEED-DEPENDENT if (iforc_tr.eq.1) then call it_catch(nt_tratm,tr_tgrid,nstep,it1,it2,clm_tscl) do i=1,npt j = ((iox(i)-1)/nx) + 1 if (yy(j) .le. -20) then ilat = 1 elseif ((yy(j) .le. 20.).and.(yy(j).gt.-20.)) then ilat = 2 elseif (yy(j) .gt. 20.) then ilat = 3 endif call val_interp(i,npt,ny,nlat_tratm,nt_tratm,ilat,it1,it2, & tr_atm,val_forw,val_back,yy) atm_val = val_back + clm_tscl*(val_forw - val_back) atm_val = 0.1*atm_val + 100.0 term_mult = 1.75/secpyr ftr(i,1,itrac) = ftr(i,1,itrac) & + term_mult*(tpw(i,1)-2)*(atm_val - tr(i,1,itrac)) enddo endif c CONSTANT WINDSPEED if (iforc_tr.eq.2) then call it_catch(nt_tratm,tr_tgrid,nstep,it1,it2,clm_tscl) do i=1,npt j = ((iox(i)-1)/nx) + 1 if (yy(j) .le. -20) then ilat = 1 elseif ((yy(j) .le. 20.).and.(yy(j).gt.-20.)) then ilat = 2 elseif (yy(j) .gt. 20.) then ilat = 3 endif call val_interp(i,npt,ny,nlat_tratm,nt_tratm,ilat,it1,it2, & tr_atm,val_forw,val_back,yy) atm_val = val_back + clm_tscl*(val_forw - val_back) atm_val = 0.1*atm_val + 100.0 ftr(i,1,itrac) = ftr(i,1,itrac) + + 50.*(atm_val -tr(i,1,itrac))*factor_c14 enddo endif c PREBOMB WITH WINDSPEED DEPENDENCE c AND W/O SUESS EFFECT c (note that setting atm_val=100. c in model units is equivalent to c setting atm_c14=0.0) if (iforc_tr.eq.3) then do i=1,npt atm_val = 100.0 term_mult = 1.75/secpyr ftr(i,1,itrac) = ftr(i,1,itrac) & + term_mult*(tpw(i,1)-2)*(atm_val - tr(i,1,itrac)) enddo endif c PREBOMB W/O WINDSPEED DEPENDENCE c AND W/O SUESS EFFECT c (note that setting atm_val=100. c in model units is equivalent to c setting atm_c14=0.0) if (iforc_tr.eq.4) then atm_val = 100. do i=1,npt ftr(i,1,itrac) = ftr(i,1,itrac) + + 50.*(atm_val-tr(i,1,itrac))*factor_c14 enddo endif if (iforc_tr.eq.5) then call it_catch(nt_tratm,tr_tgrid,nstep,it1,it2,clm_tscl) do i=1,npt j = ((iox(i)-1)/nx) + 1 if (yy(j) .le. -20) then ilat = 1 elseif ((yy(j) .le. 20.).and.(yy(j).gt.-20.)) then ilat = 2 elseif (yy(j) .gt. 20.) then ilat = 3 endif call val_interp(i,npt,ny,nlat_tratm,nt_tratm,ilat,it1,it2, & tr_atm,val_forw,val_back,yy) atm_val = val_back + clm_tscl*(val_forw - val_back) atm_val = 0.1*atm_val + 100.0 term_mult = 1.75/secpyr term_new = term_mult*(tpw(i,1)-2)*(atm_val - tr(i,1,itrac)) c SIMULATED C14 w/ windspeed dependence ftr(i,1,itrac) = ftr(i,1,itrac) + term_new c PERTURBATION C14 w/ windspeed ftr (i,1,itracp) = ftr(i,1,itracp) + term_new enddo endif if (iforc_tr.eq.6) then call it_catch(nt_tratm,tr_tgrid,nstep,it1,it2,clm_tscl) do i=1,npt j = ((iox(i)-1)/nx) + 1 if (yy(j) .le. -20) then ilat = 1 elseif ((yy(j) .le. 20.).and.(yy(j).gt.-20.)) then ilat = 2 elseif (yy(j) .gt. 20.) then ilat = 3 endif call val_interp(i,npt,ny,nlat_tratm,nt_tratm,ilat,it1,it2, & tr_atm,val_forw,val_back,yy) atm_val = val_back + clm_tscl*(val_forw - val_back) atm_val = 0.1*atm_val + 100.0 term_mult = 1.75/secpyr term_new = term_mult*(tpw(i,1)-2)*(atm_val - tr(i,1,itrac)) c SIMULATED C14 w/ windspeed dependence ftr(i,1,itrac) = ftr(i,1,itrac) + term_new c PERTURBATION C14 w/ windspeed ftr(i,1,2) = ftr(i,1,2) + term_new c AGE IN YEARS tr(i,1,3) = 0.0 ftr(i,1,3) = 0.0 do k=2,nzi(i) ftr(i,k,3) = ftr(i,k,3) + factor_age*h(i,k) enddo enddo endif if (iforc_tr.eq.7) then call it_catch(nt_tratm,tr_tgrid,nstep,it1,it2,clm_tscl) do i=1,npt j = ((iox(i)-1)/nx) + 1 if (yy(j) .le. -20) then ilat = 1 elseif ((yy(j) .le. 20.).and.(yy(j).gt.-20.)) then ilat = 2 elseif (yy(j) .gt. 20.) then ilat = 3 endif call val_interp(i,npt,ny,nlat_tratm,nt_tratm,ilat,it1,it2, & tr_atm,val_forw,val_back,yy) atm_val = val_back + clm_tscl*(val_forw - val_back) atm_val = 0.1*atm_val + 100.0 term_new = 50.*(atm_val -tr(i,1,itrac))*factor_c14 ftr(i,1,itrac) = ftr(i,1,itrac) + term_new ftr(i,1,2) = ftr(i,1,2) + term_new enddo endif if (iforc_tr.eq.8) then call it_catch(nt_tratm,tr_tgrid,nstep,it1,it2,clm_tscl) do i=1,npt j = ((iox(i)-1)/nx) + 1 if (yy(j) .le. -20) then ilat = 1 elseif ((yy(j) .le. 20.).and.(yy(j).gt.-20.)) then ilat = 2 elseif (yy(j) .gt. 20.) then ilat = 3 endif call val_interp(i,npt,ny,nlat_tratm,nt_tratm,ilat,it1,it2, & tr_atm,val_forw,val_back,yy) atm_val = val_back + clm_tscl*(val_forw - val_back) atm_val = 0.1*atm_val + 100.0 term_new = 50.*(atm_val -tr(i,1,itrac))*factor_c14 c SIMULATED C14 w/o windspeed dependence ftr(i,1,itrac) = ftr(i,1,itrac) + term_new c PERTURBATION C14 w/o windspeed dependence ftr(i,1,2) = ftr(i,1,2) + term_new c AGE IN YEARS tr(i,1,3) = 0.0 ftr(i,1,3) = 0.0 do k=2,nzi(i) ftr(i,k,3) = ftr(i,k,3) + factor_age*h(i,k) enddo enddo endif if (iforc_tr.eq.9) then call it_catch(nt_tratm,tr_tgrid,nstep,it1,it2,clm_tscl) do i=1,npt j = ((iox(i)-1)/nx) + 1 if (yy(j) .le. -20) then ilat = 1 elseif ((yy(j) .le. 20.).and.(yy(j).gt.-20.)) then ilat = 2 elseif (yy(j) .gt. 20.) then ilat = 3 endif call val_interp(i,npt,ny,nlat_tratm,nt_tratm,ilat,it1,it2, & tr_atm,val_forw,val_back,yy) atm_val = val_back + clm_tscl*(val_forw - val_back) atm_val = 0.1*atm_val + 100.0 term_mult = 1.75/secpyr term_new = term_mult*(tpw(i,1)-2)*(atm_val - tr(i,1,1)) c SIMULATED C14 w/ windspeed dependence ftr(i,1,1) = ftr(i,1,1) + term_new c PERTURBATION C14 w/ windspeed ftr(i,1,2) = ftr(i,1,2) + term_new c AGE IN YEARS tr(i,1,3) = 0.0 ftr(i,1,3) = 0.0 do k=2,nzi(i) ftr(i,k,3) = ftr(i,k,3) + factor_age*h(i,k) enddo do k=1,nzi(i) pdens_ik = pdens(i,k) c j = (iox(i)-1)/nx + 1 if (yy(j) .ge. 18.0) then if (pdens_ik.le.24.) then ftr(i,k,4)=ftr(i,k,4) & +h(i,k)*(trstar-tr(i,k,4))*factor_1month tr(i,k,5)=0.0 ftr(i,k,5)=0.0 tr(i,k,6)=0.0 ftr(i,k,6)=0.0 endif if ((pdens_ik.gt.24.).and.(pdens_ik.lt.26.5)) then ftr(i,k,5)=ftr(i,k,5) & +h(i,k)*(trstar-tr(i,k,5))*factor_1month tr(i,k,4)=0.0 ftr(i,k,4)=0.0 tr(i,k,6)=0.0 ftr(i,k,6)=0.0 endif if (pdens_ik .ge. 26.5) then ftr(i,k,6)=ftr(i,k,6) & +h(i,k)*(trstar-tr(i,k,6))*factor_1month tr(i,k,4)=0.0 ftr(i,k,4)=0.0 tr(i,k,5)=0.0 ftr(i,k,5)=0.0 endif endif if (yy(j) .le. -18.0) then if (pdens_ik.le.24.) then ftr(i,k,7)=ftr(i,k,7) & +h(i,k)*(trstar-tr(i,k,7))*factor_1month tr(i,k,8)=0.0 ftr(i,k,8)=0.0 tr(i,k,9)=0.0 ftr(i,k,9)=0.0 endif if ((pdens_ik.gt.24.).and.(pdens_ik.lt.26.5)) then ftr(i,k,8)=ftr(i,k,8) & +h(i,k)*(trstar-tr(i,k,8))*factor_1month tr(i,k,7)=0.0 ftr(i,k,7)=0.0 tr(i,k,9)=0.0 ftr(i,k,9)=0.0 endif if (pdens_ik.ge.26.5) then ftr(i,k,9)=ftr(i,k,9) & +h(i,k)*(trstar-tr(i,k,9))*factor_1month tr(i,k,7)=0.0 ftr(i,k,7)=0.0 tr(i,k,8)=0.0 ftr(i,k,8)=0.0 endif endif if (yy(j) .le. -36.0) then ftr(i,k,10)=ftr(i,k,10) & +h(i,k)*(trstar-tr(i,k,10))*factor_1month endif enddo enddo endif if (iforc_tr .eq. 51) then c AGE in YEARS c +Tracers 1 c +Tr_name ['AGE'] c +Tr_forcing 51 c +Trac_init [2] c +Trac_init_value [0.] do i=1,npt tr(i,1,itrac) = 0.0 enddo do i=1,npt do k=1,nzi(i) ftr(i,k,itrac) = ftr(i,k,itrac) + factor_age*h(i,k) enddo enddo endif if (iforc_tr .eq. 52) then c INVERSE(AGE+1) c +Tracers 1 c +Tr_name ['INV_AGE'] c +Tr_forcing 52 c +Trac_init [2] c +Trac_init_value [1.] do i=1,npt tr(i,1,itrac) = 1.0 enddo do i=1,npt do k=1,nzi(i) tr2 = tr(i,k,itrac)**2 ftr(i,k,itrac) = ftr(i,k,itrac) - tr2*factor_age*h(i,k) enddo enddo endif if (iforc_tr .eq. 53) then c RADIOACTIVE DECAY c half_life = 10 years c +Tracers 2 c +Tr_name ['TRITIUM' 'HELIUM'] c +Tr_forcing [53 0] c +Trac_init [2 2] c +Trac_init_value [1. 0.] rlambda = 10.*365.*24.*3600. do i=1,npt tr(i,1,itrac) = 1.0 tr(i,1,itracp) = 0.0 enddo decay_term = (exp(-dnt*alog(2.)/rlambda) - 1.)/dnt do i=1,npt do k=1,nzi(i) trit_sink = -h(i,k)*tr(i,k,itrac)*decay_term ftr(i,k,itrac) = ftr(i,k,itrac)-trit_sink ftr(i,k,itracp) = ftr(i,k,itracp)+trit_sink enddo enddo endif 100 continue enddo return end c------------------------------------------------------------ subroutine get_windspeed(npt,i,wind,wnsp_scalar) c------------------------------------------------------------ real wind(npt) c write(*,*) " within get_windspeed: " c write(*,*) " npt= ", npt c write(*,*) " i= ", i c write(*,*) " wind(1)= ", wind(1) wnsp_scalar = wind(i) c write(*,*) " wnsp_scalar= ", wnsp_scalar return end c------------------------------------------------------------ subroutine val_interp(i,npt,nyp,nlat_tratm,nt_tratm,ilat,it1,it2, & tr_atm,val_forw,val_back,yy) c------------------------------------------------------------ real tr_atm(nlat_tratm,nt_tratm) real yy(nyp) val_forw = tr_atm(ilat,it2) val_back = tr_atm(ilat,it1) c if (i.eq.2000) then c write(*,*) " w/i val_interp: " c write(*,*) " i= ", i c write(*,*) " val_forw= ", val_forw c write(*,*) " val_back= ", val_back c write(*,*) " npt= ", npt c write(*,*) " nyp= ", nyp c write(*,*) " nlat_tratm= ", nlat_tratm c write(*,*) " nt_tratm= ", nt_tratm c write(*,*) " ilat= ", ilat c write(*,*) " it1= ", it1 c write(*,*) " it2= ", it2 c endif return end