Given u and v wind components, compute the divergence via Spherepack.
function uv2dvF( u : float, v : float ) function uv2dvG( u : float, v : float )
Note: For the arrays whose last two dimensions are nlat x nlon, the rest of the dimensions (if any) are collectively referred to as nt. If the input/output arrays are just two dimensions, then nt can either be considered equal to 1 or nothing at all.
Arrays which have dimensions nt x nlat x nlon should not include the cyclic (wraparound) points when invoking the procedures and functions which use spherical harmonics (Spherepack).
For example, if an array x has dimensions nlat = 64 and nlon = 129, where the "129" represents the cyclic points, then the user should pass the data to the procedure/function via:
z = sample ( x([...],:,0:nlon-2) ) ; does not include cyclic points
begin nlat = 64 ; dimensions mlon = 128 mlon1 = mlon+1 fbfile = "uv300.hs" ; Generic Workstation setup nrec = fbinnumrec(fbfile) ; total number of records in the file ntim = nrec/2 ; number of time steps in dataset uvmsg = 1e+36 work = new ( (/nlat,mlon1/), float, uvmsg ) u = new ( (/nlat,mlon /), float, uvmsg ) ; source u v = new ( (/nlat,mlon /), float, uvmsg ) ; source v dvx = new ( (/nlat,mlon /), float, uvmsg ) ; divergence (same as dv) vrx = new ( (/nlat,mlon /), float, uvmsg ) ; vorticity (same as vr) sf = new ( (/nlat,mlon /), float, uvmsg ) ; stream function vp = new ( (/nlat,mlon /), float, uvmsg ) ; velocity potential uy = new ( (/nlat,mlon /), float, uvmsg ) ; latitudinal derivative (u) vy = new ( (/nlat,mlon /), float, uvmsg ) ; latitudinal derivative (v) uu = new ( (/nlat,mlon /), float, uvmsg ) ; reconstructed u vv = new ( (/nlat,mlon /), float, uvmsg ) ; reconstructed v ux = new ( (/nlat,mlon /), float, uvmsg ) ; reconstructed u vx = new ( (/nlat,mlon /), float, uvmsg ) ; reconstructed v do i = 0,nrec-1,2 month = 1 ; january if (i .ge. 2) then month = 7 ; july end if work = fbinrecread(fbfile,i ,(/nlat,mlon1/),"float") u = work(:,0:mlon-1) work = fbinrecread(fbfile,i+1,(/nlat,mlon1/),"float") v = work(:,0:mlon-1) dv = uv2dvG (u,v) ; u,v ==> divergence vr = uv2vrG (u,v) ; u,v ==> vorticity (rel) uv2vrdvg (u,v, vrx,dvx) ; u,v ==> div and vort uv2sfvpg (u,v, sf,vp) ; u,v ==> stream function + velocity pot ud = new ( (/nlat,mlon /), float, uvmsg ) vd = new ( (/nlat,mlon /), float, uvmsg ) ur = new ( (/nlat,mlon /), float, uvmsg ) dv2uvg (dv, ud,vd) ; dv ==> divergent wind components vr2uvg (vr, ur,vr) ; vr ==> rotational wind components vrdv2uvg (vr,dv, uu,vv) ; vr,dv > reconstruct original wind vrdv2uvg (vrx,dvx, ux,vx) ; vr,dv > reconstruct original wind end do end
1 : error in the specification of nlat
2 : error in the specification of nlon
4 : error in the specification of nt (jer only)
NG4.1 Home, Index, Examples, Glossary, Feedback, Ref Contents, Ref WhereAmI?