C +-======-+ C Copyright (c) 2003-2007 United States Government as represented by C the Admistrator of the National Aeronautics and Space Administration. C All Rights Reserved. C C THIS OPEN SOURCE AGREEMENT ("AGREEMENT") DEFINES THE RIGHTS OF USE, C REPRODUCTION, DISTRIBUTION, MODIFICATION AND REDISTRIBUTION OF CERTAIN C COMPUTER SOFTWARE ORIGINALLY RELEASED BY THE UNITED STATES GOVERNMENT AS C REPRESENTED BY THE GOVERNMENT AGENCY LISTED BELOW ("GOVERNMENT AGENCY"). C THE UNITED STATES GOVERNMENT, AS REPRESENTED BY GOVERNMENT AGENCY, IS AN C INTENDED THIRD-PARTY BENEFICIARY OF ALL SUBSEQUENT DISTRIBUTIONS OR C REDISTRIBUTIONS OF THE SUBJECT SOFTWARE. ANYONE WHO USES, REPRODUCES, C DISTRIBUTES, MODIFIES OR REDISTRIBUTES THE SUBJECT SOFTWARE, AS DEFINED C HEREIN, OR ANY PART THEREOF, IS, BY THAT ACTION, ACCEPTING IN FULL THE C RESPONSIBILITIES AND OBLIGATIONS CONTAINED IN THIS AGREEMENT. C C Government Agency: National Aeronautics and Space Administration C Government Agency Original Software Designation: GSC-15354-1 C Government Agency Original Software Title: GEOS-5 GCM Modeling Software C User Registration Requested. Please Visit http://opensource.gsfc.nasa.gov C Government Agency Point of Contact for Original Software: C Dale Hithon, SRA Assistant, (301) 286-2691 C C +-======-+ PROGRAM main implicit none character*120 title character*120 begdate real, allocatable :: vz(:,:) real, allocatable :: tz(:,:) real, allocatable :: wz(:,:) real, allocatable :: thz(:,:) real, allocatable :: vptpz(:,:) real, allocatable :: vpthpz(:,:) real, allocatable :: pl(:,:) real, allocatable :: pk(:,:) real, allocatable :: strm(:,:) real, allocatable :: res(:,:) real, allocatable :: vstar(:,:) real, allocatable :: wstar(:,:) real*4, allocatable :: dum(:) character*120, allocatable :: arg(:) character*120 tag, output real undef, lat0 integer im,jm,lm,tm integer i,j,L,n,nt,lrec integer iostat,rc,nargs,iargc undef = 1e15 nargs = iargc() if( nargs.ne.0 ) then allocate( arg(nargs) ) do n=1,nargs call getarg(n,arg(n)) enddo do n=1,nargs if( trim(arg(n)).eq.'-tag' ) tag = '_' // arg(n+1) enddo else tag = "" endif c Read TXT Files to determine resolution c -------------------------------------- open (10,file='LAT0' // trim(tag) // '.txt',form='formatted') read (10,*) lat0 print *, 'LAT0 = ',lat0 close(10) open (10,file='XDIM' // trim(tag) // '.txt',form='formatted') read (10,*) im print *, 'IM = ',im close(10) open (10,file='YDIM' // trim(tag) // '.txt',form='formatted') read (10,*) jm print *, 'JM = ',jm close(10) open (10,file='ZDIM' // trim(tag) // '.txt',form='formatted') read (10,*) lm print *, 'LM = ',lm close(10) open (10,file='TDIM' // trim(tag) // '.txt',form='formatted') read (10,*) tm print *, 'TDIM = ',tm close(10) open (10,file='BEGDATE' // trim(tag) // '.txt',form='formatted') read (10,*) begdate close(10) allocate( vz(jm,lm) ) allocate( tz(jm,lm) ) allocate( wz(jm,lm) ) allocate( thz(jm,lm) ) allocate( vptpz(jm,lm) ) allocate( vpthpz(jm,lm) ) allocate( pl(jm,lm) ) allocate( pk(jm,lm) ) allocate( strm(jm,lm) ) allocate( res(jm,lm) ) allocate( vstar(jm,lm) ) allocate( wstar(jm,lm) ) allocate( dum(jm) ) c GRADS Datasets c -------------- open (10,file= 'grads' // trim(tag) // '.fwrite',form='unformatted',access='direct',recl=jm*4) open (20,file='residual' // trim(tag) // '.data' ,form='unformatted',access='sequential') open (30,file='residual' // trim(tag) // '.ctl' ,form='formatted') c Read data from grads.fwrite c --------------------------- nt = 0 rc = 0 lrec = 1 do n=1,tm do L=1,lm read(10,rec=lrec,iostat=rc) dum ; vz(:,L) = dum(:) ; lrec = lrec+1 enddo do L=1,lm read(10,rec=lrec,iostat=rc) dum ; tz(:,L) = dum(:) ; lrec = lrec+1 enddo do L=1,lm read(10,rec=lrec,iostat=rc) dum ; vptpz(:,L) = dum(:) ; lrec = lrec+1 enddo do L=1,lm read(10,rec=lrec,iostat=rc) dum ; pl(:,L) = dum(:) ; lrec = lrec+1 enddo do L=1,lm read(10,rec=lrec,iostat=rc) dum ; pk(:,L) = dum(:) ; lrec = lrec+1 enddo nt = nt+1 do L=1,lm do j=1,jm if( abs(tz(j,L)-undef).gt.0.1 ) then thz(j,L) = tz(j,L)/pk(j,L) vpthpz(j,L) = vptpz(j,L)/pk(j,L) else thz(j,L) = undef vpthpz(j,L) = undef endif enddo enddo c Compute Meridional Streamfunction c --------------------------------- call stream ( vz,pl,jm,lm,strm,undef ) c Compute Mean Vertical Velocity from Continuity c ---------------------------------------------- call make_w ( vz,pl,jm,lm,wz,undef ) c Compute Residual Circulation c ---------------------------- call residual ( vz,vpthpz,thz,wz,pl,jm,lm,res,vstar,wstar,undef ) call GLAWRT (strm ,jm,LM,20) call GLAWRT (res ,jm,LM,20) call GLAWRT (vstar ,jm,LM,20) call GLAWRT (wstar ,jm,LM,20) enddo close(10) c Write Grads Control File c ------------------------ output = '^residual' // trim(tag) // '.data' title = 'Streamfunction and Residual Circulation' write(30,101) trim(output),trim(title),undef,jm,lat0,2*abs(lat0)/(jm-1),lm do L=1,lm print *, 'Pressure = ',pl(1,L) write(30,102) pl(1,L) enddo print *, 'Finished , nt = ',nt write(30,103) nt,trim(begdate),lm,lm,lm,lm 101 format('dset ',a,/, . 'title ',a,/, . 'options sequential template',/, . 'undef ',e15.6,/, . 'xdef 1 linear -180 1',/, . 'ydef ',i4,' linear ',f8.3,2x,f8.3,/, . 'zdef ',i3,' levels ') 102 format(10x,f8.3) 103 format('tdef ',i3,' linear ',a,' 1mo',/, . 'vars 4 ',/, . 'str ',i3,' 0 Streamfunction',/, . 'res ',i3,' 0 Residual Circulation',/, . 'vstar ',i3,' 0 Vstar',/, . 'wstar ',i3,' 0 wstar',/, . 'endvars') stop end SUBROUTINE ZONAL ( A,AZ,IM,JNP,undef ) DIMENSION A(IM,JNP), AZ(JNP) DO J=1,JNP AZ(J) = 0.0 IC = 0 DO I = 1, IM IF( abs(A(I,J)-UNDEF).gt.0.1 ) THEN AZ(J) = AZ(J) + A(I,J) IC = IC + 1 ENDIF ENDDO IF(IC.NE.0) AZ(J) = AZ(J) / IC IF(IC.eq.0) AZ(J) = UNDEF ENDDO RETURN END SUBROUTINE GLAWRT (A, IM,LM, KTP) real A (IM,LM) real*4 TEM (IM) DO L=1,LM DO I=1,IM if( abs(a(i,L)).gt.1.e-20 ) then TEM (I) = A(I,L) else TEM (I) = 0. endif ENDDO WRITE(KTP) TEM ENDDO RETURN END SUBROUTINE GLAWRT2 (A, JM,LM, KTP, undef) real A (JM,LM) real TEM (JM) real DPHI (JM), qout(181) real lon,lat,lons(181),lats(181) real*4 dum(181) pi = 4.0*atan(1.0) dlam = 1 dphi(:) = pi/(jm-1) iout = 1 jout = 181 dlout = 2*pi/iout dpout = pi/(jout-1) loc = 0 do j=1,jout do i=1,iout loc = loc + 1 lon = -pi + (i-1)*dlout lons(loc) = lon enddo enddo loc = 0 do j=1,jout lat = -pi/2.0 + (j-1)*dpout do i=1,iout loc = loc + 1 lats(loc) = lat enddo enddo DO L=1,LM DO J=1,JM if( abs(a(J,L)).gt.1.e-20 ) then TEM (J) = A(J,L) else TEM (J) = 0. endif ENDDO call interp_h ( tem,1,JM,1, . dlam,dphi,0.0,90.0,0.0, . qout,iout*jout,lons,lats, . 1,1,.true.,undef ) dum = qout WRITE(KTP) dum ENDDO RETURN END subroutine stream ( v0,p0,jm,lm,s,undef ) real v(jm,lm), v0(jm,lm) real s(jm,lm) real p0(jm,lm), p(jm,lm) c Invert VWND and P level index (in order to be top=>bottom) c ---------------------------------------------------------- do L=1,lm p(:,L) = p0(:,lm-L+1) v(:,L) = v0(:,lm-L+1) enddo pi = 4.*atan(1.) dp = pi/(jm-1) a = 6.372e6 g = 9.8 const = 2*pi*a/g*100 / 1.e10 do k=1,lm do j=1,jm s(j,k) = 0.0 enddo do L=1,k do j=1,jm delp = (p(j,L+1)-p(j,L-1))/2 if(L.eq.1) delp = (p(j,L+1)-p(j,L) )/2 if(L.eq.lm) delp = (p(j,L) -p(j,L-1))/2 phi = -pi/2+(j-1)*dp if( abs(v(j,L)-undef).gt.0.1 ) then s(j,k) = s(j,k) + v(j,l)*delp*cos(phi)*const endif enddo enddo enddo c Invert Streamfunction for grads output (in order to be bottom=>top) c ------------------------------------------------------------------- do k=1,lm do j=1,jm v(j,k) = s(j,lm-k+1) enddo enddo do k=1,lm do j=1,jm s(j,k) = v(j,k) enddo enddo return end subroutine residual ( v0,vpthp0,th0,w0,p0,jm,lm,res,vstar,wstar,undef ) real v0(jm,lm), v(jm,lm) real w0(jm,lm), w(jm,lm) real vpthp0(jm,lm), th0(jm,lm) real vpthp (jm,lm), th (jm,lm), dthdp(jm,lm) real stuff (jm,lm) real res (jm,lm) real vstar(jm,lm) real wstar(jm,lm) real s (jm,lm) real p0(jm,lm), p(jm,lm), rho0(jm,lm) real cosp(jm) real ddcosp(jm,lm) logical defined PARAMETER ( AIRMW = 28.97 ) PARAMETER ( RUNIV = 8314.3 ) PARAMETER ( CPD = 1004.16 ) PARAMETER ( RGAS = RUNIV/AIRMW) PARAMETER ( AKAP = RGAS/CPD ) c Invert v,th,vpthp, and P level index (in order to be top=>bottom) c ----------------------------------------------------------------- do L=1,lm w(:,L) = w0(:,lm-L+1) v(:,L) = v0(:,lm-L+1) p(:,L) = p0(:,lm-L+1) th(:,L) = th0(:,lm-L+1) vpthp(:,L) = vpthp0(:,lm-L+1) enddo pi = 4.*atan(1.) dp = pi/(jm-1) a = 6.372e6 g = 9.8 H = 7000.0 ps = 1000.0 ts = 240.0 rhos = ps/(rgas*ts) c Compute Mean Air Density c ------------------------ do L=1,lm do j=1,jm z = -H*log(p(j,L)/ps) rho0(j,L) = rhos*exp(-z/H) enddo enddo c Compute D(Theta)/DZ (with a forced minimum to prevent dthdz => 0) c ----------------------------------------------------------------- do L=2,lm-1 do j=1,jm delp = p(j,L+1)-p(j,L-1) if( defined(th(j,L+1),undef ) .and. . defined(th(j,L-1),undef ) ) then dthdp(j,L) = min( -0.003, ( th(j,L+1)-th(j,L-1) )/ delp ) else dthdp(j,L) = undef endif enddo enddo do j=1,jm dthdp(j,1) = dthdp(j,2) dthdp(j,lm) = dthdp(j,lm-1) enddo c Compute D(cos*vpthp/dthdz)/Dphi c ------------------------------- do L=1,lm do j=1,jm phi = -pi/2 + (j-1)*dp cosp(j) = cos(phi) if( defined(vpthp(j,L),undef) .and. . defined(dthdp(j,L),undef) .and. . dthdp(j,L).ne.0.0 ) then stuff(j,L) = -H*cosp(j)*vpthp(j,L)/(p(j,L)*dthdp(j,L)) else stuff(j,L) = undef endif enddo enddo do L=1,lm do j=1,jm if( j.gt.1 .and. j.lt.jm ) then if( defined(stuff(j+1,L),undef) .and. . defined(stuff(j-1,L),undef) ) then ddcosp(j,L) = ( stuff(j+1,L)-stuff(j-1,L) )/(2*dp) else ddcosp(j,L) = undef endif else ddcosp(j,L) = undef endif enddo enddo c Compute Wstar c ------------- do L=1,lm do j=1,jm if( defined(ddcosp(j,L),undef) ) then wstar(j,Lm-L+1) = w(j,L) + ddcosp(j,L)/(a*cosp(j)) else wstar(j,Lm-L+1) = undef endif enddo enddo c Compute D(rho*vpthp/dthdz)/DZ c ----------------------------- do L=1,lm do j=1,jm if( defined(dthdp(j,L),undef) .and. . defined(vpthp(j,L),undef) .and. . dthdp(j,L).ne.0.0 ) then stuff(j,L) = rho0(j,L)*vpthp(j,L)/(p(j,L)*dthdp(j,L)) else stuff(j,L) = undef endif enddo enddo do L=2,lm-1 do j=1,jm delp = p(j,L+1)-p(j,L-1) if( defined(stuff(j,L+1),undef) .and. . defined(stuff(j,L-1),undef) ) then res(j,L) = ( stuff(j,L+1)-stuff(j,L-1) )/ delp else res(j,L) = undef endif enddo enddo do j=1,jm res(j,1) = res(j,2) res(j,lm) = res(j,lm-1) enddo c Compute Vstar c ------------- do L=1,lm do j=1,jm if( defined( res(j,L),undef) .and. . defined( v(j,L),undef) ) then res(j,L) = v(j,L) - p(j,L)*res(j,L)/rho0(j,L) else res(j,L) = undef endif vstar(j,LM-L+1) = res(j,L) enddo enddo c Construct Streamfunction for Vstar c ---------------------------------- do k=1,lm do j=1,jm s(j,k) = 0.0 enddo do L=1,k do j=1,jm delp = (p(j,L+1)-p(j,L-1))/2 if(L.eq.1) delp = (p(j,L+1)-p(j,L) )/2 if(L.eq.lm) delp = (p(j,L) -p(j,L-1))/2 if( defined(res(j,L),undef) ) then s(j,k) = s(j,k) + res(j,L)*delp*cosp(j)*rho0(j,L)*H/p(j,L) endif enddo enddo enddo c Invert Streamfunction for grads output (in order to be bottom=>top) c ------------------------------------------------------------------- do L=1,lm do j=1,jm res(j,L) = s(j,LM-L+1) enddo enddo return end subroutine make_w ( v0,p0,jm,lm,w,undef ) real v0(jm,lm), v(jm,lm) real p0(jm,lm), p(jm,lm) real w(jm,lm), rho0(jm,lm) real s(jm,lm), cosp(jm) real dvcos_dphi(jm,lm) logical defined PARAMETER ( AIRMW = 28.97 ) PARAMETER ( RUNIV = 8314.3 ) PARAMETER ( CPD = 1004.16 ) PARAMETER ( RGAS = RUNIV/AIRMW) PARAMETER ( AKAP = RGAS/CPD ) c Invert v and P level index (in order to be top=>bottom) c ------------------------------------------------------- do L=1,lm v(:,L) = v0(:,lm-L+1) p(:,L) = p0(:,lm-L+1) enddo pi = 4.*atan(1.) dp = pi/(jm-1) a = 6.372e6 g = 9.8 H = 7000.0 ps = 1000.0 ts = 240.0 rhos = ps/(rgas*ts) do L=1,lm do j=1,jm phi = -pi/2 + (j-1)*dp cosp(j) = cos(phi) z = -H*log(p(j,L)/ps) rho0(j,L) = rhos*exp(-z/H) enddo enddo do L=1,lm do j=1,jm if( j.gt.1 .and. j.lt.jm ) then if( defined(v(j+1,L),undef) .and. . defined(v(j-1,L),undef) ) then dvcos_dphi(j,L) = ( v(j+1,L)*cosp(j+1)-v(j-1,L)*cosp(j-1) )/(2*dp) else dvcos_dphi(j,L) = undef endif else dvcos_dphi(j,L) = undef endif enddo enddo c Construct W from Continuity c --------------------------- do k=1,lm do j=1,jm s(j,k) = undef enddo do L=1,k do j=1,jm delp = (p(j,L+1)-p(j,L-1))/2 if(L.eq.1) delp = (p(j,L+1)-p(j,L) )/2 if(L.eq.lm) delp = (p(j,L) -p(j,L-1))/2 phi = -pi/2+(j-1)*dp if( dvcos_dphi(j,L).ne.undef ) then if( s(j,k).eq.undef ) s(j,k) = 0.0 s(j,k) = s(j,k) + dvcos_dphi(j,L)*delp*rho0(j,L)*H/(p(j,L)*a*cosp(j)) endif enddo enddo do j=1,jm if( s(j,k).ne.undef ) s(j,k) = s(j,k)/rho0(j,k) enddo enddo c Invert Streamfunction for grads output (in order to be bottom=>top) c ------------------------------------------------------------------- do k=1,lm do j=1,jm w(j,k) = s(j,lm-k+1) enddo enddo return end subroutine interp_h ( q_cmp,im,jm,lm, . dlam,dphi,rotation,tilt,precession, . q_geo,irun,lon_geo,lat_geo, . msgn,norder,check,undef ) C*********************************************************************** C C PURPOSE: C ======== C Performs a horizontal interpolation from a field on a computational grid C to arbitrary locations. C C INPUT: C ====== C q_cmp ...... Field q_cmp(im,jm,lm) on the computational grid C im ......... Longitudinal dimension of q_cmp C jm ......... Latitudinal dimension of q_cmp C lm ......... Vertical dimension of q_cmp C dlam ....... Computational Grid Delta Lambda C dphi ....... Computational Grid Delta Phi C rotation ... Rotation parameter lam_np (Degrees) C tilt ....... Rotation parameter phi_np (Degrees) C precession . Rotation parameter lam_0 (Degrees) C irun ....... Number of Output Locations C lon_geo .... Longitude Location of Output C lat_geo .... Latitude Location of Output C msgn ....... Flag for scalar field ( msgn = 1 ) C or vector component ( msgn = -1 ) C norder ..... Order of Interpolation: Bi-Linear => abs(norder) = 1 C Bi-Cubic => abs(norder) = 3 C Note: If norder < 0, then check for positive definite C check ...... Logical Flag to check for Undefined values C C OUTPUT: C ======= C q_geo ...... Field q_geo(irun,lm) at arbitrary locations C C C*********************************************************************** C* GODDARD LABORATORY FOR ATMOSPHERES * C*********************************************************************** implicit none c Input Variables c --------------- integer im,jm,lm,irun,norder,msgn logical check real q_geo(irun,lm) real lon_geo(irun) real lat_geo(irun) real q_cmp(im,jm,lm) real dlam(im) real dphi(jm) c Local Variables c --------------- integer i,j,l,m,n integer, allocatable :: ip0(:), im1(:) integer, allocatable :: jp0(:), jm1(:) c Bi-Linear Weights c ----------------- real, allocatable :: wl_ip0jp0 (:) real, allocatable :: wl_im1jp0 (:) real, allocatable :: wl_ip0jm1 (:) real, allocatable :: wl_im1jm1 (:) real, allocatable :: old_lon (:) real, allocatable :: old_lat (:) real, allocatable :: old_dlam(:) real, allocatable :: old_dphi(:) real ap0, am1 real bp0, bm1 real lon_cmp(im) real lat_cmp(jm) real q_tmp(irun) real pi,cosnp,sinnp,p1,p2,p3,eps,d real lam,lam_ip1,lam_ip0,lam_im1,lam_im2 real phi,phi_jp1,phi_jp0,phi_jm1,phi_jm2 real dl,dp,lam_np,phi_np,lam_0,eps_np real rotation , tilt , precession real lam_geo, lam_cmp real phi_geo, phi_cmp real undef, getcon integer im1_cmp,icmp integer jm1_cmp,jcmp logical compute_weights real old_rotation real old_tilt real old_precession data old_rotation /-999.9/ data old_tilt /-999.9/ data old_precession /-999.9/ parameter ( eps = 1.e-10 ) c Initialization c -------------- pi = 4.*atan(1.) dl = 2*pi/ im ! Uniform Grid Delta Lambda dp = pi/(jm-1) ! Uniform Grid Delta Phi c Allocate Memory for Weights and Index Locations c ----------------------------------------------- if(.not.allocated(old_lon)) then allocate ( old_dlam(im) , old_dphi(jm) ) allocate ( old_lon(irun) , old_lat(irun) ) allocate ( wl_ip0jp0(irun) , wl_im1jp0(irun) ) allocate ( wl_ip0jm1(irun) , wl_im1jm1(irun) ) allocate ( ip0(irun) , im1(irun) ) allocate ( jp0(irun) , jm1(irun) ) do i=1,irun old_lon(i) = -999.9 old_lat(i) = -999.9 enddo do i=1,im old_dlam(i) = 0.0 enddo do j=1,jm old_dphi(j) = 0.0 enddo else i = size (old_dlam) j = size (old_dphi) m = size (old_lon) if(i.ne.im .or. j.ne.jm .or. m.ne.irun) then deallocate ( old_dlam , old_dphi ) deallocate ( old_lon , old_lat ) deallocate ( wl_ip0jp0 , wl_im1jp0 ) deallocate ( wl_ip0jm1 , wl_im1jm1 ) deallocate ( ip0 , im1 ) deallocate ( jp0 , jm1 ) allocate ( old_dlam(im) , old_dphi(jm) ) allocate ( old_lon(irun) , old_lat(irun) ) allocate ( wl_ip0jp0(irun) , wl_im1jp0(irun) ) allocate ( wl_ip0jm1(irun) , wl_im1jm1(irun) ) allocate ( ip0(irun) , im1(irun) ) allocate ( jp0(irun) , jm1(irun) ) do i=1,irun old_lon(i) = -999.9 old_lat(i) = -999.9 enddo do i=1,im old_dlam(i) = 0.0 enddo do j=1,jm old_dphi(j) = 0.0 enddo endif endif c Compute Input Computational-Grid Latitude and Longitude Locations c ----------------------------------------------------------------- lon_cmp(1) = -pi do i=2,im lon_cmp(i) = lon_cmp(i-1) + dlam(i-1) enddo lat_cmp(1) = -pi*0.5 do j=2,jm-1 lat_cmp(j) = lat_cmp(j-1) + dphi(j-1) enddo lat_cmp(jm) = pi*0.5 c Check for Co-incident Grid-Point Latitude and Pole Locations c ------------------------------------------------------------ eps_np = 0.0 do j=1,jm phi_cmp = lat_cmp(j)*180./pi if( abs( phi_cmp-tilt ).lt.1.e-3 ) eps_np = 1.e-3 if( tilt+eps_np .gt. 90. ) eps_np = -1.e-3 enddo lam_np = pi/180.*rotation phi_np = pi/180.*(tilt+eps_np) lam_0 = pi/180.*precession if( tilt.eq.90. ) then cosnp = 0.0 sinnp = 1.0 else if(tilt.eq.-90.0) then cosnp = 0.0 sinnp =-1.0 else cosnp = cos(phi_np) sinnp = sin(phi_np) endif c Determine if Weights Need to be Updated c --------------------------------------- compute_weights = rotation.ne.old_rotation .or. . tilt.ne.old_tilt .or. . precession.ne.old_precession m = 1 do while ( .not.compute_weights .and. m.le.irun ) compute_weights = (lon_geo(m).ne.old_lon(m)) .or. . (lat_geo(m).ne.old_lat(m)) m = m+1 enddo i = 1 do while ( .not.compute_weights .and. i.le.im ) compute_weights = dlam(i).ne.old_dlam(i) i = i+1 enddo j = 1 do while ( .not.compute_weights .and. j.le.jm-1 ) compute_weights = dphi(j).ne.old_dphi(j) j = j+1 enddo c Compute Weights for Computational to Geophysical Grid Interpolation c ------------------------------------------------------------------- if( compute_weights ) then old_rotation = rotation old_tilt = tilt old_precession = precession do i=1,irun old_lon(i) = lon_geo(i) old_lat(i) = lat_geo(i) lam_geo = lon_geo(i) phi_geo = lat_geo(i) p1 = cosnp*cos(phi_geo)*cos(lam_geo+lam_0-pi) . + sin(phi_geo)*sinnp p1 = min(p1, 1.0) p1 = max(p1,-1.0) phi_cmp = asin( p1 ) if( tilt.eq.90.0 .or. tilt.eq.-90.0 ) then p2 = sinnp*cos(lam_geo+lam_0-pi) else p2 = sinnp*cos(phi_geo)*cos(lam_geo+lam_0-pi) . - sin(phi_geo)*cosnp p2 = p2 / max( cos(phi_cmp),eps ) p2 = min(p2, 1.0) p2 = max(p2,-1.0) endif p2 = acos( p2 ) p3 = cos(phi_geo)*sin(lam_geo+lam_0-pi) if( p3.lt.0.0 ) p2 = -p2 p2 = p2 + lam_np - pi lam_cmp = mod( p2+3.0*pi,2.0*pi ) - pi c Determine Indexing Based on Computational Grid c ---------------------------------------------- im1_cmp = 1 do icmp = 2,im if( lon_cmp(icmp).lt.lam_cmp ) im1_cmp = icmp enddo jm1_cmp = 1 do jcmp = 2,jm if( lat_cmp(jcmp).lt.phi_cmp ) jm1_cmp = jcmp enddo im1(i) = im1_cmp ip0(i) = im1(i) + 1 jm1(i) = jm1_cmp jp0(i) = jm1(i) + 1 c Fix Longitude Index Boundaries c ------------------------------ if(im1(i).eq.im) then ip0(i) = 1 endif c Compute Immediate Surrounding Coordinates c ----------------------------------------- lam = lam_cmp phi = phi_cmp c Compute and Adjust Longitude Weights c ------------------------------------ lam_im1 = lon_cmp(im1(i)) lam_ip0 = lon_cmp(ip0(i)) if( lam_im1.gt.lam_ip0 ) lam_ip0 = lam_ip0 + 2*pi c Compute and Adjust Latitude Weights c Note: Latitude Index Boundaries are Adjusted during Interpolation c ------------------------------------------------------------------ phi_jm1 = lat_cmp(jm1(i)) phi_jp0 = lat_cmp(jp0(i)) if( jm1(i).eq.jm ) phi_jp0 = phi_jm1 + dphi(jm-1) c Linear Weights c -------------- d = (phi_jp0-phi_jm1) wl_ip0jm1(i) = (phi_jp0-phi )/d wl_ip0jp0(i) = (phi -phi_jm1)/d enddo endif c Interpolate Computational-Grid Quantities to Geophysical Grid Using Bi-Linear c ----------------------------------------------------------------------------- do L=1,lm do i=1,irun if( q_cmp( ip0(i),jm1(i),L ).ne.undef .and. . q_cmp( ip0(i),jp0(i),L ).ne.undef ) then q_tmp(i) = wl_ip0jm1(i) * q_cmp( ip0(i),jm1(i),L ) . + wl_ip0jp0(i) * q_cmp( ip0(i),jp0(i),L ) else q_tmp(i) = undef endif enddo c Load Temp array into Output array c --------------------------------- do i=1,irun q_geo(i,L) = q_tmp(i) enddo enddo deallocate ( old_dlam , old_dphi ) deallocate ( old_lon , old_lat ) deallocate ( wl_ip0jp0 , wl_im1jp0 ) deallocate ( wl_ip0jm1 , wl_im1jm1 ) deallocate ( ip0 , im1 ) deallocate ( jp0 , jm1 ) return end function defined ( q,undef ) implicit none logical defined real q,undef defined = abs(q-undef).gt.0.1*undef return end function defined