c c------------------------------------------------------------------ c Begin file /home/sting/foster/loc/mkwi.f c------------------------------------------------------------------ c subroutine mkwi(plotp,idim1,idim2,loc,iht,iparm) c c if iparm = 1 -> return unmag (do not calc ion drifts) c if iparm = 2 -> return vnmag (do not calc ion drifts) c if iparm = 3 -> Calculate ui, return in plotp c if iparm = 4 -> Calculate vi, return in plotp c if iparm = 5 -> Calculate wi, return in plotp c If iht <= 0 -> find iparm at pressure surfaces c If iht > 0 -> find iparm at constant height surfaces c include 'tgcmparam.h' include 'input.h' include 'tigcmfld.h' include 'loc.h' dimension plotp(idim1,idim2) dimension tn(kmx),un(kmx),vn(kmx),xo1(kmx),xo2(kmx),xn2(kmx), + z(kmx),wn(kmx) c c if iefield > 0 -> call efield in wi calculation data iseasav/0/, iutav/0/ save iseasav,iutav c if (iparm.lt.1.or.iparm.gt.nextra) then write(6,"('mkwi: bad iparm=',i5,' mkwi returning')") iparm return endif dipr = dipdec(1,loc)*dtr decr = dipdec(2,loc)*dtr si = sin(dipr) sii = si*si ci = cos(dipr) cii = ci*ci bgaus = 0.4 xmaso = 32. woi = 9.57946e+3 * bgaus / xmaso colfac = 1.5 cosd = cos(decr) sind = sin(decr) day = iyd(1) - iyd(1)/1000*1000 xlag=0. c c Find wi at gcm pressure grid surfaces (use array utzp for gcm fields): if (iht.le.0) then do 100 it=1,ntms c c Call efield, if necessary (all velocities in m/s): if (iefield.gt.0.and.iparm.le.5.and.iparm.ge.3) then c call efield(gmloc(1,loc),gmloc(2,loc),day,ut(it),iseasav, c + iutav,pot,vu,ve) timeut=ut(it)+xlag call efield(gmloc(1,loc),gmloc(2,loc),day,timeut, + iseasav,iutav,pot,vu,ve) uit = ve vut = vu else uit = 0. vut = 0. endif do 105 k=1,kmx c vinp = 0.73e-9 * sqrt(utzp(it,k,ixt,loc)/1000.) * c + utzp(it,k,ixo1,loc) * colfac + c + 0.69e-9 * utzp(it,k,ixn2,loc) + c + 0.67e-9 * utzp(it,k,ixo2,loc) vinp=3.e-10*(utzp(it,k,ixo1,loc)+utzp(it,k,ixn2,loc) + +utzp(it,k,ixo2,loc)) rhoi = vinp/woi coi = 1. / (1. + rhoi*rhoi) unmag = utzp(it,k,ixu,loc)*cosd - utzp(it,k,ixv,loc)*sind vnmag = utzp(it,k,ixu,loc)*sind + utzp(it,k,ixv,loc)*cosd wn1 = utzp(it,k,ixw,loc) c unmag=0. c vnmag=0. c wn1=0. if (iparm.eq.1) then ! unmag plotp(it,k) = unmag elseif (iparm.eq.2) then ! vnmag plotp(it,k) = vnmag elseif (iparm.eq.3) then ! ui plotp(it,k)=coi*(rhoi*rhoi*unmag+rhoi*(-vnmag*si-wn1*ci)) + +coi*(uit+rhoi*vut) elseif (iparm.eq.4) then ! vi plotp(it,k)=coi*(rhoi*rhoi*vnmag+rhoi*unmag*si+ + vnmag*cii-wn1*si*ci)+coi*(vut*si-rhoi*uit*si) elseif (iparm.eq.5) then ! wi plotp(it,k)=coi*(rhoi*rhoi*wn1+rhoi*unmag*ci-vnmag*si*ci+ + wn1*sii)+coi*(vut*ci-rhoi*uit*ci) c plotp(it,k) = (-vnmag*ci*si + rhoi*unmag*ci) * coi endif 105 continue 100 continue c c Find wi at nhtscale(nhtscale) constant height surfaces: c (interpolate utzp gcm fields) c (return cpspval in plotp if height out of range) c else do 200 it=1,ntms tn(:) = utzp(it,:,ixt,loc) un(:) = utzp(it,:,ixu,loc) vn(:) = utzp(it,:,ixv,loc) wn(:) = utzp(it,:,ixw,loc) xo1(:) = utzp(it,:,ixo1,loc) xo2(:) = utzp(it,:,ixo2,loc) xn2(:) = utzp(it,:,ixn2,loc) z(:) = utzp(it,:,ixz,loc) if (iefield.gt.0.and.iparm.le.5.and.iparm.ge.3) then c call efield(gmloc(1,loc),gmloc(2,loc),day,ut(it),iseasav, c + iutav,pot,vu,ve) timeut=ut(it)+xlag call efield(gmloc(1,loc),gmloc(2,loc),day,timeut, + iseasav,iutav,pot,vu,ve) uit = ve vut = vu else uit = 0. vut = 0. endif do 205 k=1,nhtscale tnht = terpo(tn,z,kmx,htscale(k),0,cpspval) unht = terpo(un,z,kmx,htscale(k),0,cpspval) vnht = terpo(vn,z,kmx,htscale(k),0,cpspval) wnht = terpo(wn,z,kmx,htscale(k),0,cpspval) xo1ht = terpo(xo1,z,kmx,htscale(k),1,cpspval) xo2ht = terpo(xo2,z,kmx,htscale(k),1,cpspval) xn2ht = terpo(xn2,z,kmx,htscale(k),1,cpspval) if (tnht.eq.cpspval.or.unht.eq.cpspval.or.vnht.eq.cpspval + .or.wnht.eq.cpspval.or.xo1ht.eq.cpspval + .or.xo2ht.eq.cpspval.or.xn2ht.eq.cpspval) then plotp(it,k) = cpspval else c vinp = 0.73e-9 * sqrt(tnht/1000.) * xo1ht * colfac + c + 0.69e-9 * xn2ht + 0.67e-9 * xo2ht vinp=3.e-10*(utzp(it,k,ixo1,loc)+utzp(it,k,ixn2,loc) + +utzp(it,k,ixo2,loc)) rhoi = vinp/woi coi = 1. / (1. + rhoi*rhoi) unmag = unht*cosd - vnht*sind vnmag = unht*sind + vnht*cosd c unmag=0. c vnmag=0. c wnht=0. if (iparm.eq.1) then ! unmag plotp(it,k) = unmag elseif (iparm.eq.2) then ! vnmag plotp(it,k) = vnmag elseif (iparm.eq.3) then ! ui plotp(it,k)=coi*(rhoi*rhoi*unmag+rhoi* + (-vnmag*si-wnht*ci))+coi*(uit+rhoi*vut) elseif (iparm.eq.4) then ! vi plotp(it,k)=coi*(rhoi*rhoi*vnmag+rhoi*unmag*si+ + vnmag*cii-wnht*si*ci)+coi*(vut*si- + rhoi*uit*si) elseif (iparm.eq.5) then ! wi plotp(it,k)=coi*(rhoi*rhoi*wnht+rhoi*unmag*ci- + vnmag*si*ci+wnht*sii)+coi*(vut*ci-rhoi*uit*ci) c plotp(it,k) = (-vnmag*ci*si + rhoi*unmag*ci) * coi endif endif 205 continue 200 continue endif c return end c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c function terpo(var,varht,idim,ht,log,spval) c c Interpolate var to height ht, using column of heights in varht c if log <= 0 -> linear interpolation c if log > 0 -> exponential interpolation c (Return spval if ht not in varht) c dimension var(idim),varht(idim) c c Bracket the height: do 100 i=2,idim if (ht.ge.varht(i-1).and.ht.le.varht(i)) then i1 = i-1 i2 = i goto 101 endif 100 continue c write(6,"('terpo: could not find ht=',f10.3,' in array', c + ' varht=',/(8f10.3))") ht,varht c write(6,"('terpo returning spval=',e12.4)") spval terpo = spval return 101 continue c c Linear interpolation: if (log.le.0) then frac = (ht-varht(i1)) / (varht(i2)-varht(i1)) terpo = frac*var(i2) + (1.-frac)*var(i1) c c Exponential (log) interpolation: else terpo = var(i1) * exp((alog(var(i2) / var(i1)) / + (varht(i2)-varht(i1))) * (ht-varht(i1))) endif c return end