!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& Program metrcs ! ****************************************************************** ! * * ! * Convert Map Coordinates into Cartesian * ! * * ! ****************************************************************** implicit none !---------------------------------------------------------------------- include 'param_o.h' integer,Parameter::i0=im+im-1 integer,Parameter::i1=i0+i0-1 integer,Parameter::img=i1+i0-1 & , j0=i0,j1=i1,jmg=img & , iml=i0,jml=iml,lmg=img & , n1grid=img/2,nagrid=n1grid/3,ntay=30,nn=128 real,parameter::flon0=270. & , flat0=-90.,scen=1 real::JUMBLE(nn),WORK(nn),RA(0:nn),QA(0:nn) integer,dimension(jmg)::imin,imax real,dimension(lmg)::s real,dimension(3)::xe real,dimension(0:img+1,0:img+1,2)::xg,yg,zg real,dimension(0:iml+1,0:iml+1,nm)::xa,ya,za real,dimension(0:im+1,0:im+1,nm)::xh,yh,zh,xv,yv,zv integer,dimension(0:im+1,0:jm+1,nm)::imod logical,dimension(0:im+1,0:jm+1,nm)::lv real,Dimension(0:im+1,0:jm+1,nm)::DXVDX,DXVDY,DYVDX,DYVDY, & DZVDX,DZVDY real::qs11,qs12,qs22,sqhc,dxhdx,dyhdx,dzhdx,dxhdy,dyhdy,dzhdy real,dimension(0:im+1,0:jm+1,nm)::sqv,sqh,rsqv,rsqh,q11,q12,q22 real,dimension(0:im+1,0:jm+1,nm,4)::qd11,qd12,qd21,qd22 real,dimension(0:im+1,0:jm+1,nm)::qbv11,qbv12,qbv22 real,dimension(0:im+1,0:jm+1,nm)::vcosp,vsinp,vlam,vphi, & hcosp,hsinp,hlam,hphi !----------------------------------------------------------------------- !----------------------------------------------------------------------- real,Dimension(igm,jgm):: GLON , GLAT, HG, & XMP0, YMP0 integer,dimension(igm,jgm)::KMP0 !----------------------------------------------------------------------- real,Dimension(0:im+1, 0:jm+1,nm):: u2us,v2us , & u2vs, v2vs,qh11,qh12,qh22 !----------------------------------------------------------------------- Logical NCOR (0:im + 1, 0:jm + 1) real,dimension(0:im+1,0:jm+1,nm,4,2,2)::qvh,qhv,qvh2,qhv2 !----------------------------------------------------------------------- Character * 01 cres Character * 37 cfile1 Character * 37 cfile2 real:: mfac,dmap,xmap,ymap,atang1,csgm,sngm real::pihlf,pi,rad,dlat,dlon,xlon,ylat,csal,snal,csbt & ,snbt,x0,y0,z0,x1,y1,z1,x2,y2,z2,vlamda,vcosl,vsinl,ax,ay,bx & ,by,det,rdet integer::j,i,l,il,jl,n,imd,k,it,jt,nres,list1,list2,kmap CHARACTER (len=39) :: tfile=& '../data_in/grid/SOCT0004.DAT' real,dimension(2)::xm real,dimension(0:im+1,0:jm+1,nm)::rarea real,dimension(0:im+1,0:jm+1,nm,4)::ds real,parameter::a=6.371e6 !----------------------------------------------------------------------- mfac=0.5 !----------------------------------------------------------------------- ! ! Initialize grid generation ! CALL insoct2(flon0, flat0, scen, tfile) ! Call Inoct(n1grid,nagrid,ntay,flon0,flat0,scen,RA,QA,nn & ! ,JUMBLE,WORK) ! ! Define control indices for global arrays ! Do j=1,j0-1 IMIN(j)=i0+1-j IMAX(j)=i1+j-1 End Do Do j=j0,j1 IMIN(j)=1 IMAX(j)=img End Do Do j=j1+1,jmg IMIN(j)=j-j1+1 IMAX(j)=img-j+j1 End Do xm(1)=-0.13 xm(2)=-0.87 call Otoc1(xm,xe) call Otoc2(xm,xe) call ctoo(xe,xm,kmap) ! ! Read_in Global Absolute Coordinates (XG,YG,ZG) ! dmap=2./(img-1) S(1)=-1 Do l=2,lmg S(l)=S(l-1)+dmap End Do Do j=1,jmg Do i=IMIN(j),IMAX(j) xm(1)=S(i) xm(2)=S(j) Call Otoc1(xm,XE) XG(i,j,1)=XE(1) YG(i,j,1)=XE(2) ZG(i,j,1)=XE(3) Call Otoc2(xm,XE) XG(i,j,2)=XE(1) YG(i,j,2)=XE(2) ZG(i,j,2)=XE(3) End Do End Do ! ! Define Local Arrays in the Corners of the Northern Hemisphere ! ! Domain 6 ! Do j=1,j0 Do i=IMIN(j),i0 Xa(i,j,6)=XG(i,j,1) Ya(i,j,6)=YG(i,j,1) Za(i,j,6)=ZG(i,j,1) End Do End Do Do j=1,j0-1 Do i=1,IMIN(j)-1 Xa(i,j,6)=XG(i1+j-1,j0+1-i,2) Ya(i,j,6)=YG(i1+j-1,j0+1-i,2) Za(i,j,6)=ZG(i1+j-1,j0+1-i,2) End Do End Do ! ! Domain 8 ! Do j=1,j0 Do i=i1,IMAX(j) il=i-i1+1 Xa(il,j,8)=XG(i,j,1) Ya(il,j,8)=YG(i,j,1) Za(il,j,8)=ZG(i,j,1) End Do End Do Do j=1,j0-1 Do i=IMAX(j)+1,img il=i-i1+1 Xa(il,j,8)=XG(i0+1-j,i-i1+1,2) Ya(il,j,8)=YG(i0+1-j,i-i1+1,2) Za(il,j,8)=ZG(i0+1-j,i-i1+1,2) End Do End Do ! ! Domain 12 ! Do j=j1,jmg Do i=IMIN(j),i0 jl=j-j1+1 Xa(i,jl,12)=XG(i,j,1) Ya(i,jl,12)=YG(i,j,1) Za(i,jl,12)=ZG(i,j,1) End Do End Do Do j=j1+1,jmg Do i=1,IMIN(j)-1 jl=j-j1+1 Xa(i,jl,12)=XG(i1+img-j,j1+i-1,2) Ya(i,jl,12)=YG(i1+img-j,j1+i-1,2) Za(i,jl,12)=ZG(i1+img-j,j1+i-1,2) End Do End Do ! ! Domain 14 ! Do j=j1,jmg Do i=i1,IMAX(j) il=i-i1+1 jl=j-j1+1 Xa(il,jl,14)=XG(i,j,1) Ya(il,jl,14)=YG(i,j,1) Za(il,jl,14)=ZG(i,j,1) End Do End Do Do j=j1+1,jmg Do i=IMAX(j)+1,img il=i-i1+1 jl=j-j1+1 Xa(il,jl,14)=XG(jl,jmg-il+1,2) Ya(il,jl,14)=YG(jl,jmg-il+1,2) Za(il,jl,14)=ZG(jl,jmg-il+1,2) End Do End Do ! ! Define All Other Local Arrays ! Xa(1:iml,1:jml,1)=XG(i0:i1,1:j0,2) Ya(1:iml,1:jml,1)=YG(i0:i1,1:j0,2) Za(1:iml,1:jml,1)=ZG(i0:i1,1:j0,2) Xa(1:iml,1:jml,2)=XG(1:i0,j0:j1,2) Ya(1:iml,1:jml,2)=YG(1:i0,j0:j1,2) Za(1:iml,1:jml,2)=ZG(1:i0,j0:j1,2) Xa(1:iml,1:jml,3)=XG(i0:i1,j0:j1,2) Ya(1:iml,1:jml,3)=YG(i0:i1,j0:j1,2) Za(1:iml,1:jml,3)=ZG(i0:i1,j0:j1,2) Xa(1:iml,1:jml,4)=XG(i1:img,j0:j1,2) Ya(1:iml,1:jml,4)=YG(i1:img,j0:j1,2) Za(1:iml,1:jml,4)=ZG(i1:img,j0:j1,2) Xa(1:iml,1:jml,5)=XG(i0:i1,j1:jmg,2) Ya(1:iml,1:jml,5)=YG(i0:i1,j1:jmg,2) Za(1:iml,1:jml,5)=ZG(i0:i1,j1:jmg,2) Xa(1:iml,1:jml,7)=XG(i0:i1,1:j0,1) Ya(1:iml,1:jml,7)=YG(i0:i1,1:j0,1) Za(1:iml,1:jml,7)=ZG(i0:i1,1:j0,1) Xa(1:iml,1:jml,9)=XG(1:i0,j0:j1,1) Ya(1:iml,1:jml,9)=YG(1:i0,j0:j1,1) Za(1:iml,1:jml,9)=ZG(1:i0,j0:j1,1) Xa(1:iml,1:jml,10)=XG(i0:i1,j0:j1,1) Ya(1:iml,1:jml,10)=YG(i0:i1,j0:j1,1) Za(1:iml,1:jml,10)=ZG(i0:i1,j0:j1,1) Xa(1:iml,1:jml,11)=XG(i1:img,j0:j1,1) Ya(1:iml,1:jml,11)=YG(i1:img,j0:j1,1) Za(1:iml,1:jml,11)=ZG(i1:img,j0:j1,1) Xa(1:iml,1:jml,13)=XG(i0:i1,j1:jmg,1) Ya(1:iml,1:jml,13)=YG(i0:i1,j1:jmg,1) Za(1:iml,1:jml,13)=ZG(i0:i1,j1:jmg,1) Call Bocoh(Xa,iml,jml) Call Bocoh(Ya,iml,jml) Call Bocoh(Za,iml,jml) csal=Cos(alfa) snal=Sin(alfa) csbt=Cos(beta) snbt=Sin(beta) csgm=cos(gamm) sngm=sin(gamm) do n=1,nm do i=0,iml+1 do j=0,iml+1 x0=xa(i,j,n) y0=ya(i,j,n) z0=za(i,j,n) x1= x0*csal+y0*snal y1=-x0*snal+y0*csal z1=z0 x2=x1 y2= y1*csbt+z1*snbt z2=-y1*snbt+z1*csbt x1= x2*csgm+y2*sngm y1=-x2*sngm+y2*csgm z1=z2 xa(i,j,n)=x1 ya(i,j,n)=y1 za(i,j,n)=z1 enddo enddo enddo ! ! Define XH,YH,ZH ! Do n=1,14 XH(1:im,1:jm,n)=Xa(1:iml:2,1:jml:2,n) YH(1:im,1:jm,n)=Ya(1:iml:2,1:jml:2,n) ZH(1:im,1:jm,n)=Za(1:iml:2,1:jml:2,n) End Do Call Bocoh(XH,im,jm) Call Bocoh(YH,im,jm) Call Bocoh(ZH,im,jm) do n=1,nm do i=0,im+1 do j=0,jm+1 HCOSP (i, j, n) = Sqrt (Xh (i,j, n) **2 + Yh (i,j, n) **2) HSINP (i, j, n) = Zh (i,j, n) HLAM (i, j, n) = Atang1 (Yh (i,j, n), Xh (i,j, n) ) HPHI (i, j, n) = ASIN (HSINP (i, j, n) ) enddo enddo enddo ! print *,'center of tile 3,',hlam((im-1)/2,(jm-1)/2,3)*180./3.14, & ! hphi((im-1)/2,(jm-1)/2,3)*180./3.14,gamm ! stop ! print *,hlam(9,13,9),yh(9,13,9),xh(9,13,9) ! stop ! ! Define XV,YV,ZV ! Do n=1,nm XV(0:im,0:jm,n)=Xa(0:iml+1:2,0:jml+1:2,n) YV(0:im,0:jm,n)=Ya(0:iml+1:2,0:jml+1:2,n) ZV(0:im,0:jm,n)=Za(0:iml+1:2,0:jml+1:2,n) End Do do n=1,nm do i=0,im do j=0,jm VCOSP (i, j, n) = Sqrt (Xv (i,j, n) **2 + Yv (i,j, n) **2) VSINP (i, j, n) = Zv (i,j, n) VLAM (i, j, n) = Atang1 (Yv (i,j, n), Xv (i,j, n) ) VPHI (i, j, n) = ASIN (VSINP (i, j, n) ) enddo enddo enddo ! ! Correction of coefficients in the corners ! IMOD(:,:,:)=1 IMOD(1,1, 1)=0 IMOD(1,1, 2)=0 IMOD(1,1, 7)=0 IMOD(1,1, 8)=0 IMOD(1,1, 9)=0 IMOD(1,1,12)=0 IMOD(im,1,6 )=0 IMOD(im,1,7 )=0 IMOD(im,1,11)=0 IMOD(im,1,14)=0 IMOD(im,1, 1)=0 IMOD(1,jm, 2)=0 IMOD(im,1, 4)=0 IMOD(1,jm, 5)=0 IMOD(1,jm, 6)=0 IMOD(1,jm, 9)=0 IMOD(1,jm,13)=0 IMOD(1,jm,14)=0 IMOD(im,jm, 4)=0 IMOD(im,jm, 5)=0 IMOD(im,jm, 8)=0 IMOD(im,jm,11)=0 IMOD(im,jm,12)=0 IMOD(im,jm,13)=0 LV(:,:,:)=.true. LV(0 ,0 ,1)=.false. LV(im,0 ,1)=.false. LV(0 ,0 ,2)=.false. LV(0 ,jm,2)=.false. LV(0 ,jm,5)=.false. LV(im,jm,5)=.false. LV(im,jm,4)=.false. LV(im,0 ,4)=.false. LV(0 ,0 ,8)=.false. LV(im,jm,8)=.false. LV(im,0 ,14)=.false. LV(0 ,jm,14)=.false. LV(im,jm,12)=.false. LV(0 ,0 ,12)=.false. LV(im,0 ,6 )=.false. LV(0 ,jm,6 )=.false. LV(0 ,0 ,7)=.false. LV(im,0 ,7)=.false. LV(0 ,0 ,9)=.false. LV(0 ,jm,9)=.false. LV(0 ,jm,13)=.false. LV(im,jm,13)=.false. LV(im,jm,11)=.false. LV(im,0 ,11)=.false. Do n = 1, nm Do j = 0, jm Do i = 0, im DXVDX (i,j,n) =(XH(i+1,j,n)-XH(i,j+1,n))*mfac DYVDX (i,j,n) =(YH(i+1,j,n)-YH(i,j+1,n))*mfac DZVDX (i,j,n) =(ZH(i+1,j,n)-ZH(i,j+1,n))*mfac DXVDY (i,j,n) =(XH(i+1,j+1,n)-XH(i,j,n))*mfac DYVDY (i,j,n) =(YH(i+1,j+1,n)-YH(i,j,n))*mfac DZVDY (i,j,n) =(ZH(i+1,j+1,n)-ZH(i,j,n))*mfac qs11 = DXVDX (i,j,n)**2+DYVDX(i,j,n)**2+DZVDX(i,j,n)**2 qs22 = DXVDY(i,j,n)**2+DYVDY(i,j,n)**2+DZVDY(i,j,n)**2 qs12 = DXVDX (i, j, n) * DXVDY (i, j, n) + DYVDX (i, j, n) & * DYVDY (i, j, n) + DZVDX (i, j, n) * DZVDY (i, j, n) If (LV(i, j,n) ) Then SQV (i, j,n) = Sqrt (qs11 * qs22 - qs12 * qs12) RSQV (i, j,n) = 1. / SQV (i, j,n) Q11 (i, j, n) = qs22 * RSQV (i, j,n) **2 Q22 (i, j, n) = qs11 * RSQV (i, j,n) **2 Q12 (i, j, n) = - qs12 * RSQV (i, j,n) **2 else SQV (i, j,n) = 0. RSQV (i, j,n) = 0. Q11 (i, j, n) = 0. Q22 (i, j, n) = 0. Q12 (i, j, n) = 0. EndIf ! print *,i,j,n,sqv(i,j,n),qs11,qs22,qs12 EndDo EndDo EndDo call bocov1(sqv,im,jm) call bocov2(q11,q22,im,jm) call bocov3(q12,im,jm) call qvhqhv(xh,yh,zh,xv,yv,zv,q11,q12,q22,qvh,qhv,im,jm,nm,mfac) call qvhqhv2(xh,yh,zh,xv,yv,zv,q11,q12,q22,qh11,qh12,qh22,qvh2 & ,qhv2,im,jm,nm,mfac) call qintc(xh,yh,zh,xv,yv,zv,q11,q12,q22,im,jm,nm,mfac) !DRAGAN call qintc_2(xh,yh,zh,xv,yv,zv,im,jm,nm,mfac) ! ! Grid components: H-points ! sqhc=SQV(1,1,7)*2./3. do n=1,nm Do j = 1, jm Do i = 1, im imd=IMOD(i,j,n) dxhdx=(XV(i,j-1,n)-XV(i-1,j,n))*mfac dyhdx=(YV(i,j-1,n)-YV(i-1,j,n))*mfac dzhdx=(ZV(i,j-1,n)-ZV(i-1,j,n))*mfac dxhdy=(XV(i,j,n)-XV(i-1,j-1,n))*mfac dyhdy=(YV(i,j,n)-YV(i-1,j-1,n))*mfac dzhdy=(ZV(i,j,n)-ZV(i-1,j-1,n))*mfac qs11 = dxhdx**2 + dyhdx**2 + dzhdx**2 qs22 = dxhdy**2 + dyhdy**2 + dzhdy**2 qs12 = dxhdx * dxhdy + dyhdx * dyhdy + dzhdx * dzhdy ! SQH(i,j,n)=Sqrt(qs11*qs22-qs12*qs12)*imd & ! + sqhc*(1-imd) sqh(i,j,n)=0.25*(sqv(i-1,j-1,n)+sqv(i,j,n)+sqv(i,j-1,n) & +sqv(i-1,j,n)) ! print *,i,j,n,sqh(i,j,n) rarea(i,j,n)=1./(2.*sqh(i,j,n)*a*a) ds(i,j,n,1)=0.5*sqrt(qs11)*a ds(i,j,n,3)=0.5*sqrt(qs22)*a EndDo EndDo EndDo call bocoh(sqh,im,jm) call bocoh(rarea,im,jm) do n=1,nm Do j = 1, jm Do i = 1, im dxhdx=(XV(i,j-1,n)+XV(i,j,n)-XV(i-1,j,n)-XV(i-1,j-1,n))*0.5 dyhdx=(YV(i,j-1,n)+YV(i,j,n)-YV(i-1,j,n)-YV(i-1,j-1,n))*0.5 dzhdx=(ZV(i,j-1,n)+ZV(i,j,n)-ZV(i-1,j,n)-ZV(i-1,j-1,n))*0.5 dxhdy=(XV(i,j,n)+XV(i-1,j,n)-XV(i,j-1,n)-XV(i-1,j-1,n))*0.5 dyhdy=(YV(i,j,n)+YV(i-1,j,n)-YV(i,j-1,n)-YV(i-1,j-1,n))*0.5 dzhdy=(ZV(i,j,n)+ZV(i-1,j,n)-ZV(i,j-1,n)-ZV(i-1,j-1,n))*0.5 qs11 = dxhdx**2 + dyhdx**2 + dzhdx**2 qs22 = dxhdy**2 + dyhdy**2 + dzhdy**2 ds(i,j,n,2)=0.25*sqrt(qs11)*a ds(i,j,n,4)=0.25*sqrt(qs22)*a EndDo EndDo EndDo do n=1,nm do i=1,im1 do j=1,jm1 do k=1,4 if(k.eq.1)then it=i-1 jt=j else if(k.eq.2)then it=i+1 jt=j else if(k.eq.3)then it=i jt=j-1 else it=i jt=j+1 endif qd11(i,j,n,k)=q11(it,jt,n)*(dxvdx(it,jt,n)*dxvdx(i,j,n) & +dyvdx(it,jt,n)*dyvdx(i,j,n)+dzvdx(it,jt,n)*dzvdx(i,j,n)) & +q12(it,jt,n)*(dxvdy(it,jt,n)*dxvdx(i,j,n) & +dyvdy(it,jt,n)*dyvdx(i,j,n)+dzvdy(it,jt,n)*dzvdx(i,j,n)) qd12(i,j,n,k)=q12(it,jt,n)*(dxvdx(it,jt,n)*dxvdx(i,j,n) & +dyvdx(it,jt,n)*dyvdx(i,j,n)+dzvdx(it,jt,n)*dzvdx(i,j,n)) & +q22(it,jt,n)*(dxvdy(it,jt,n)*dxvdx(i,j,n) & +dyvdy(it,jt,n)*dyvdx(i,j,n)+dzvdy(it,jt,n)*dzvdx(i,j,n)) qd21(i,j,n,k)=q11(it,jt,n)*(dxvdx(it,jt,n)*dxvdy(i,j,n) & +dyvdx(it,jt,n)*dyvdy(i,j,n)+dzvdx(it,jt,n)*dzvdy(i,j,n)) & +q12(it,jt,n)*(dxvdy(it,jt,n)*dxvdy(i,j,n) & +dyvdy(it,jt,n)*dyvdy(i,j,n)+dzvdy(it,jt,n)*dzvdy(i,j,n)) qd22(i,j,n,k)=q12(it,jt,n)*(dxvdx(it,jt,n)*dxvdy(i,j,n) & +dyvdx(it,jt,n)*dyvdy(i,j,n)+dzvdx(it,jt,n)*dzvdy(i,j,n)) & +q22(it,jt,n)*(dxvdy(it,jt,n)*dxvdy(i,j,n) & +dyvdy(it,jt,n)*dyvdy(i,j,n)+dzvdy(it,jt,n)*dzvdy(i,j,n)) enddo enddo enddo enddo do n=1,nm do i=0,im+1 do j=0,jm+1 if(sqh(i,j,n).ne.0)then rsqh(i,j,n)=1./sqh(i,j,n) else rsqh(i,j,n)=0. endif enddo enddo enddo do n=1,nm do i=0,im do j=0,jm qbv11(i,j,n)=0.5*(q11(i,j,n)+q22(i,j,n)+2.*q12(i,j,n)) qbv12(i,j,n)=0.5*(q22(i,j,n)-q11(i,j,n)) qbv22(i,j,n)=0.5*(q11(i,j,n)+q22(i,j,n)-2.*q12(i,j,n)) ! print *,i,j,n,qbv11(i,j,n),qbv12(i,j,n),qbv22(i,j,n) enddo enddo enddo ! ! Outfile for Absolute Coordintes ! print *,"hsinp(1,1,1),",hsinp(1,1,1) list1 = 11 nres=1 write (cres, 900) nres 900 format(i1.1) cfile1 = "../data_in/grid/gridinit.dat" open (unit = list1, file = cfile1, form = "unformatted") Write (list1) SQV, SQH, Q11, Q12, Q22,DXVDX, DXVDY, DYVDX, & DYVDY, DZVDX, DZVDY, RSQV, VCOSP, VSINP, VLAM, VPHI, RSQH, HCOSP, & HSINP, HLAM, HPHI,qd11,qd12,qd21,qd22,qvh,qhv,qbv11,qbv12,qbv22 & ,qvh2,qhv2,qh11,qh12,qh22 Close (list1) open(list1,file="../data_in/grid/hvphi.dat", & form="unformatted") write(list1)hphi,vphi close(list1) open(list1,file="../data_in/grid/hpos.dat", & form="unformatted") write(list1)hlam,hphi close(list1) open(list1,file="../data_in/grid/fzeff.dat", & form="unformatted") write(list1)rarea,ds close(list1) !test print * , 'Wrote ', cfile1, ' file' print *,'center of tile 3,',hlam((im-1)/2,(jm-1)/2,3), & hphi((im-1)/2,(jm-1)/2,3) !test !----------------------------------------------------------------------- ! ! Prepare Graphical Output ! ! ! Define GLON and GLAT ! pihlf=2.*Atan(1.) pi=2.*pihlf rad=pi/180. dlat= pi/(jgm-1) dlon=2 *pi/igm Do i=1,igm xlon=(i-1)*dlon Do j=1,jgm GLON(i,j)=xlon End Do End Do Do j=1,jgm ylat = - pihlf + (j - 1) * dlat Do i=1,igm GLAT(i,j)=ylat End Do End Do csal=Cos(alfa) snal=Sin(alfa) csbt=Cos(beta) snbt=Sin(beta) ! ! Find XMP0, YMP0, KMP0 ! Do j=1,jgm Do i=1,igm xlon=GLON(i,j) ylat=GLAT(i,j) x2=Cos(ylat)*Cos(xlon) y2=Cos(ylat)*Sin(xlon) z2=Sin(ylat) x0 = x2 * csgm - y2 * sngm y0 = x2 * sngm + y2 * csgm z0 = z2 x1=x0 y1=y0*csbt-z0*snbt z1=y0*snbt+z0*csbt x0=x1*csal-y1*snal y0=x1*snal+y1*csal z0=z1 XE(1)=x0 XE(2)=y0 XE(3)=z0 Call Ctoo(XE,xm,kmap) xmap=xm(1) ymap=xm(2) xmap=(1.+xmap)/2. ymap=(1.+ymap)/2. XMP0(i,j)=xmap YMP0(i,j)=ymap KMP0(i,j)=kmap ! print *,i,j,xmap,ymap,kmap ! stop End Do End Do ! ! Prepare conversion of covariant back to spherical winds ! do n = 1, nm do j = 1, jm1 do i = 1, im1 vlamda = VLAM (i, j, n) vcosl = cos (vlamda) vsinl = sin (vlamda) ax = - vsinl * DXVDX (i, j, n) + vcosl * DYVDX (i, j, n) ay = - vsinl * DXVDY (i, j, n) + vcosl * DYVDY (i, j, n) bx = - vcosl * VSINP (i, j, n) * DXVDX (i, j, n) - vsinl * VSINP ( & i, j, n) * DYVDX (i, j, n) + VCOSP (i, j, n) * DZVDX (i, j, n) by = - vcosl * VSINP (i, j, n) * DXVDY (i, j, n) - vsinl * VSINP ( & i, j, n) * DYVDY (i, j, n) + VCOSP (i, j, n) * DZVDY (i, j, n) det = ax * by - ay * bx rdet = 1. / det u2us (i, j, n) = by * rdet v2us (i, j, n) = - bx * rdet u2vs (i, j, n) = - ay * rdet v2vs (i, j, n) = ax * rdet !print *,i,j,n,u2us(i,j,n),v2us(i,j,n),u2vs(i,j,n),v2vs(i,j,n) enddo enddo enddo !----------------------------------------------------------------------- ! ! Outfile for graphics ! print *,igm,jgm,im,jm,nm list2 = 12 cfile2 = "../data_in/grid/grph.dat" open (unit = list2, file = cfile2, form = "unformatted") Write (list2) GLON, GLAT, XMP0, YMP0, KMP0, u2us, v2us, u2vs, & v2vs Close (list2) !test print * , 'Wrote ', cfile2, ' file' ! open(11,file="ocdt.dat",form="unformatted") ! write(11)xa ! write(11)ya ! write(11)za ! close(11) !test !----------------------------------------------------------------------- Stop END program metrcs