c subroutine proclat(j,iden,gcmzp,ipflds) c c Process a latitude slice for timegcm c Raw lat slice at j from history is defined in f on input c processed fields are returned in f c ipflds(nflds+2) (29+n2+fof2): if ipflds(ip) > 0 -> process field ip c include 'gettimes.h' c parameter(spv=1.e36) parameter(p0=5.e-4, boltz=1.38e-16) dimension work(imxp1,kmx),pzps(kmx),gcmzp(kmx),ipflds(nflds+2) c c The fields are as follows (n2 not on history): c (ch4 and o21d added 1/92): c tn un vn o2 o n4s noz co co2 h2o h2 hox o+ ch4 o21d c 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 c no2 no o3 npo1 oh ho2 h n2d ti te ne o2+ w z n2 c 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 c c Fields 1-25,27 are defined at 1/2 levels, so body (k=2->kmx-1) c may be linearly interpolated (thats all fields except ne,w,z): c itxn2 = nflds+1 do 100 ip=1,27 if (iden.gt.0) then if ((ipflds(itxox).gt.0.and.(ip.eq.itxo1.or.ip.eq.itxo3)).or. + (ipflds(itxnoz).gt.0.and.(ip.eq.itxno.or.ip.eq.itxno2))) + goto 101 if (ipflds(itxw).gt.0.and.ip.eq.itxz) goto 101 if (ip.eq.itxt.or.ip.eq.itxo1.or.ip.eq.itxo2.or.ip.eq.itxn2) + goto 101 endif if (ipflds(ip).le.0.or.ip.eq.itxne) goto 100 101 continue do k=2,kmx-1 f(:,k,ip) = 0.5*(f(:,k,ip)+f(:,k-1,ip)) enddo 100 continue c c Fields t,u,v have bottom level in top slot (so 1 <- kmx). c Then value at kmx gets value at kmx-1: c do ip=1,3 if (ipflds(ip).gt.0) then f(:,1,ip) = f(:,kmx,ip) f(:,kmx,ip) = f(:,kmx-1,ip) endif enddo c c Convert u,v to m/s: c do ip=2,3 if (ipflds(ip).gt.0) + f(:,:,ip) = f(:,:,ip) * .01 enddo c c Fields o2 and ox get linear extrap at bottom, linear interp at top: c do ip=4,5 if (ipflds(ip).gt.0) then f(:,1,ip) = 1.5*f(:,1,ip)-0.5*f(:,2,ip) f(:,kmx,ip) = 0.5*(f(:,kmx,ip)+f(:,kmx-1,ip)) endif enddo c c N2 is calculated as 1.-o2-ox: c gcmn2(:,:) = amax1(.00001,(1.-f(:,:,itxo2)-f(:,:,itxo1))) c c Prognostics n4s,noz,co,co2,h2o,h2,hox, and ch4 level 1 gets bottom c half level (i.e., when kmx=37, -13 <- -12.75). Then use linear c interpolation at the top: c tn un vn o2 o n4s noz co co2 h2o h2 hox o+ ch4 o21d c 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 c no2 no o3 npo1 oh ho2 h n2d ti te ne o2+ w z n2 c 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 c do ip=6,14 if (ipflds(ip).le.0.or.ip.eq.itxop) goto 150 f(:,kmx,ip) = 0.5*(f(:,kmx-1,ip)+f(:,kmx,ip)) 150 continue enddo c c Diagnostics o+,o21d,no2,no,o3,o1,oh,ho2,h,n2d,ti,te, and o2p get c bottom half level (as above), and log extrapolation at top: c do 200 ip=13,27 if (iden.gt.0) then if ((ipflds(itxox).gt.0.and.(ip.eq.itxo1.or.ip.eq.itxo3)).or. + (ipflds(itxnoz).gt.0.and.(ip.eq.itxno.or.ip.eq.itxno2))) + goto 200 if (ip.eq.itxt.or.ip.eq.itxo1.or.ip.eq.itxo2) goto 200 endif if (ip.eq.itxch4.or.ip.eq.itxne) goto 200 f(:,kmx,ip) = sqrt(f(:,kmx-1,ip)**3 / f(:,kmx-2,ip)) 200 continue c c Fields w, z, and ne are defined at full levels: c (convert z to km and add 62 km) c Note, however, that glbmean histories zpht are in km (not cm), so c z will be wrong when this code reads a source history made c by the mksrc code. (the timegcm histories are in cm, so we c leave the below conversion intact, and ignore z for source c history reads) c 9/19/91: no longer add 62 km to heights c 2/92: Changed mksrc to write cm, so this km conversion should be c ok for source histories as well as model c c Calculate scale heights and vertical velocity: c if (ipflds(itxw).gt.0) then do i=1,imx+1 do k=2,kmx-1 pzps(k) = f(i,k+1,itxz) - f(i,k-1,itxz) enddo pzps(1) = 4.*f(i,2,itxz) - 3.*f(i,1,itxz) - f(i,3,itxz) pzps(kmx) =3.*f(i,kmx,itxz) - 4.*f(i,kmx-1,itxz)+ + f(i,kmx-2,itxz) f(i,:,itxw) = (f(i,:,itxw) * pzps(:)) enddo endif if (ipflds(itxz).gt.0) f(:,:,itxz) = f(:,:,itxz) * 1.e-5 if (ipflds(itxw).gt.0) f(:,:,itxw) = f(:,:,itxw) * .01 c c Convert to number densities, if desired: c if (iden.eq.1.or.iden.eq.2) then c c For ox, work = rmox = (o1+o3) / (o1/16. + o3/48.) c if (ipflds(itxox).gt.0) then work(:,:) = (f(:,:,itxo1) + f(:,:,itxo3)) / + (f(:,:,itxo1)/gcmwtime(itxo1) + + f(:,:,itxo3) / gcmwtime(itxo3)) f(:,:,itxox) = f(:,:,itxox) * gcmwtime(itxox) / work(:,:) endif c c TEMP fix for zero bottom boundaries of no,no2 (4/8/92): c eps = 1.e-20 if (ipflds(itxnoz).gt.0) then do i=1,imxp1 do k=1,kmx if (f(i,k,itxno).lt.eps.and.f(i,k,itxno2).lt.eps) then f(i,k,itxnoz) = spv goto 98 endif work(i,k) = (f(i,k,itxno) + f(i,k,itxno2)) / + (f(i,k,itxno)/gcmwtime(itxno) + + f(i,k,itxno2) / gcmwtime(itxno2)) f(i,k,itxnoz) = f(i,k,itxnoz) * gcmwtime(itxnoz) / work(i,k) 98 continue enddo enddo endif c c TEMP fix for zero bottom boundaries of oh,ho2,h (4/8/92): c if (ipflds(itxhox).gt.0) then do i=1,imxp1 do k=1,kmx if (f(i,k,itxoh).lt.eps.and.f(i,k,itxho2).lt.eps.and. + f(i,k,itxh).lt.eps) then f(i,k,itxhox) = spv goto 99 endif work(i,k) = (f(i,k,itxoh)+f(i,k,itxho2)+ + f(i,k,itxh)) / (f(i,k,itxoh) / + gcmwtime(itxoh) + f(i,k,itxho2)/ + gcmwtime(itxho2) + + f(i,k,itxh)/gcmwtime(itxh)) f(i,k,itxhox) = f(i,k,itxhox) * gcmwtime(itxhox) / work(i,k) 99 continue enddo enddo endif c c work = p0 * e(-z) * barm / kT c do k=1,kmx work(:,k) = p0 * exp(-gcmzp(k)) / + ((f(:,k,itxo2)/gcmwtime(itxo2)+f(:,k,itxo1)/ + gcmwtime(itxo1)+gcmn2(:,k)/28.) * boltz * f(:,k,itxt)) enddo c c Now do actual conversion to number densities: c do ip=1,nflds if ((ipflds(ip).gt.0.and.gcmwtime(ip).gt.0.).or. + ip.eq.itxo1.or.ip.eq.itxo2.or.ip.eq.itxn2) + f(:,:,ip) = f(:,:,ip) * work(:,:) / gcmwtime(ip) enddo if (ipflds(itxn2).gt.0) gcmn2(:,:) = gcmn2(:,:) * work(:,:) / 28. c c Convert number densities to mixing ratios, if desired: c endif if (iden.eq.2) then work(:,:) = f(:,:,itxo2) + f(:,:,itxo1) + + gcmn2(:,:) do ip=1,nflds if (ipflds(ip).gt.0.and.gcmwtime(ip).gt.0) + f(:,:,ip) = f(:,:,ip) / work(:,:) enddo if (ipflds(itxn2).gt.0) gcmn2(:,:) = gcmn2(:,:) / work(:,:) endif c return end