grid_tools  1.6.0
 All Data Structures Files Functions Variables
pesg.f90
Go to the documentation of this file.
1 
4 
15 module pesg
16 !=============================================================================
17 use pkind, only: spi,dp
18 use pietc, only: f,t,u0,u1,u2,o2,rtod,dtor,pih,pi2
19 implicit none
20 private
21 public :: xctoxm_ak,xmtoxc_ak,get_edges,bestesg_geo,bestesg_map, &
22  hgrid_ak_rr,hgrid_ak_rc,hgrid_ak_dd,hgrid_ak_dc,hgrid_ak, &
23  gtoxm_ak_rr,gtoxm_ak_dd,xmtog_ak_rr,xmtog_ak_dd
24 
25 interface xctoxs; module procedure xctoxs; end interface
26 interface xstoxc; module procedure xstoxc,xstoxc1; end interface
27 interface xstoxt; module procedure xstoxt; end interface
28 interface xttoxs; module procedure xttoxs,xttoxs1; end interface
29 interface xttoxm; module procedure xttoxm; end interface
30 interface zttozm; module procedure zttozm; end interface
31 interface xmtoxt; module procedure xmtoxt,xmtoxt1; end interface
32 interface zmtozt; module procedure zmtozt,zmtozt1; end interface
33 interface xctoxm_ak; module procedure xctoxm_ak; end interface
34 interface xmtoxc_ak
35  module procedure xmtoxc_ak,xmtoxc_vak,xmtoxc_vak1; end interface
36 interface get_edges; module procedure get_edges; end interface
37 interface get_qx; module procedure get_qx,get_qxd; end interface
38 interface get_qofv;module procedure get_qofv,get_qofvd,get_qsofvs;end interface
39 interface get_meanq; module procedure get_meanqd,get_meanqs; end interface
40 interface guessak_map; module procedure guessak_map; end interface
41 interface guessak_geo; module procedure guessak_geo; end interface
42 interface bestesg_geo; module procedure bestesg_geo; end interface
43 interface bestesg_map; module procedure bestesg_map; end interface
44 interface hgrid_ak_rr;module procedure hgrid_ak_rr,hgrid_ak_rr_c; end interface
45 interface hgrid_ak_rc; module procedure hgrid_ak_rc; end interface
46 interface hgrid_ak_dd;module procedure hgrid_ak_dd,hgrid_ak_dd_c; end interface
47 interface hgrid_ak_dc; module procedure hgrid_ak_dc; end interface
48 interface hgrid_ak; module procedure hgrid_ak,hgrid_ak_c; end interface
49 interface gtoxm_ak_rr
50  module procedure gtoxm_ak_rr_m,gtoxm_ak_rr_g; end interface
51 interface gtoxm_ak_dd
52  module procedure gtoxm_ak_dd_m,gtoxm_ak_dd_g; end interface
53 interface xmtog_ak_rr
54  module procedure xmtog_ak_rr_m,xmtog_ak_rr_g; end interface
55 interface xmtog_ak_dd
56  module procedure xmtog_ak_dd_m,xmtog_ak_dd_g; end interface
57 
58 interface gaulegh; module procedure gaulegh; end interface
59 
60 contains
61 
67 subroutine xctoxs(xc,xs)! [xctoxs]
68 implicit none
69 real(dp),dimension(3),intent(in ):: xc
70 real(dp),dimension(2),intent(out):: xs
71 real(dp):: zp
72 zp=u1+xc(3); xs=xc(1:2)/zp
73 end subroutine xctoxs
74 
84 subroutine xstoxc(xs,xc,xcd)! [xstoxc]
85 use pmat4, only: outer_product
86 implicit none
87 real(dp),dimension(2), intent(in ):: xs
88 real(dp),dimension(3), intent(out):: xc
89 real(dp),dimension(3,2),intent(out):: xcd
90 real(dp):: zp
91 zp=u2/(u1+dot_product(xs,xs)); xc(1:2)=xs*zp; xc(3)=zp
92 xcd=-outer_product(xc,xs)*zp; xcd(1,1)=xcd(1,1)+zp; xcd(2,2)=xcd(2,2)+zp
93 xc(3)=xc(3)-u1
94 end subroutine xstoxc
95 
106 subroutine xstoxc1(xs,xc,xcd,xcdd)! [xstoxc]
107 use pmat4, only: outer_product
108 implicit none
109 real(dp),dimension(2), intent(in ):: xs
110 real(dp),dimension(3), intent(out):: xc
111 real(dp),dimension(3,2), intent(out):: xcd
112 real(dp),dimension(3,2,2),intent(out):: xcdd
113 real(dp),dimension(3,2):: zpxcdxs
114 real(dp),dimension(3) :: zpxc
115 real(dp) :: zp
116 integer(spi) :: i
117 zp=u2/(u1+dot_product(xs,xs)); xc(1:2)=xs*zp; xc(3)=zp
118 xcd=-outer_product(xc,xs)*zp
119 zpxc=zp*xc; xc(3)=xc(3)-u1; xcdd=u0
120 do i=1,2
121  zpxcdxs=xcd*xc(i)
122  xcdd(:,i,i)=xcdd(:,i,i)-zpxc
123  xcdd(:,i,:)=xcdd(:,i,:)-zpxcdxs
124  xcdd(:,:,i)=xcdd(:,:,i)-zpxcdxs
125  xcdd(i,:,i)=xcdd(i,:,i)-zpxc(1:2)
126  xcdd(i,i,:)=xcdd(i,i,:)-zpxc(1:2)
127 enddo
128 do i=1,2; xcd(i,i)=xcd(i,i)+zp; enddo
129 end subroutine xstoxc1
130 
138 subroutine xstoxt(k,xs,xt,ff)! [xstoxt]
139 implicit none
140 real(dp), intent(in ):: k
141 real(dp),dimension(2),intent(in ):: xs
142 real(dp),dimension(2),intent(out):: xt
143 logical, intent(out):: ff
144 real(dp):: s,sc
145 s=k*(xs(1)*xs(1)+xs(2)*xs(2)); sc=u1-s
146 ff=abs(s)>=u1; if(ff)return
147 xt=u2*xs/sc
148 end subroutine xstoxt
149 
158 subroutine xttoxs(k,xt,xs,xsd,ff)! [xttoxs
159 use pmat4, only: outer_product
160 implicit none
161 real(dp), intent(in ):: k
162 real(dp),dimension(2), intent(in ):: xt
163 real(dp),dimension(2), intent(out):: xs
164 real(dp),dimension(2,2),intent(out):: xsd
165 logical, intent(out):: ff
166 real(dp),dimension(2):: rspd
167 real(dp) :: s,sp,rsp,rsppi,rsppis
168 integer(spi) :: i
169 s=k*dot_product(xt,xt); sp=u1+s
170 ff=(sp<=u0); if(ff)return
171 rsp=sqrt(sp)
172 rsppi=u1/(u1+rsp)
173 rsppis=rsppi**2
174 xs=xt*rsppi
175 rspd=k*xt/rsp
176 xsd=-outer_product(xt,rspd)*rsppis
177 do i=1,2; xsd(i,i)=xsd(i,i)+rsppi; enddo
178 end subroutine xttoxs
179 
192 subroutine xttoxs1(k,xt,xs,xsd,xsdd,xs1,xsd1,ff)! [xttoxs]
193 use pmat4, only: outer_product
194 implicit none
195 real(dp), intent(in ):: k
196 real(dp),dimension(2), intent(in ):: xt
197 real(dp),dimension(2), intent(out):: xs ,xs1
198 real(dp),dimension(2,2), intent(out):: xsd,xsd1
199 real(dp),dimension(2,2,2),intent(out):: xsdd
200 logical, intent(out):: ff
201 real(dp),dimension(2,2):: rspdd
202 real(dp),dimension(2) :: rspd,rspd1,rsppid
203 real(dp) :: s,sp,rsp,rsppi,rsppis,s1,rsp1
204 integer(spi) :: i
205 s1=dot_product(xt,xt); s=k*s1; sp=u1+s
206 ff=(sp<=u0); if(ff)return
207 rsp=sqrt(sp); rsp1=o2*s1/rsp
208 rsppi=u1/(u1+rsp); rsppis=rsppi**2
209 xs=xt*rsppi; xs1=-xt*rsp1*rsppis
210 rspd=k*xt/rsp; rspd1=(xt*rsp-k*xt*rsp1)/sp
211 rsppid=-rspd*rsppis
212 xsd1=-outer_product(xt,rspd1-u2*rspd*rsp1*rsppi)
213 do i=1,2; xsd1(i,i)=xsd1(i,i)-rsp1; enddo; xsd1=xsd1*rsppis
214 
215 xsd=-outer_product(xt,rspd)*rsppis
216 do i=1,2; xsd(i,i)=xsd(i,i)+rsppi; enddo
217 
218 rspdd=-outer_product(xt,rspd)*rsppi
219 xsdd=u0
220 do i=1,2; xsdd(i,:,i)= rsppid; enddo
221 do i=1,2; xsdd(i,i,:)=xsdd(i,i,:)+rsppid; enddo
222 do i=1,2; xsdd(:,:,i)=xsdd(:,:,i)+u2*rspdd*rsppid(i); enddo
223 do i=1,2; rspdd(i,i)=rspdd(i,i)+rsp*rsppi; enddo
224 do i=1,2; xsdd(i,:,:)=xsdd(i,:,:)-xt(i)*rspdd*rsppi*k/sp; enddo
225 end subroutine xttoxs1
226 
234 subroutine xttoxm(a,xt,xm,ff)! [xttoxm]
235 implicit none
236 real(dp), intent(in ):: a
237 real(dp),dimension(2),intent(in ):: xt
238 real(dp),dimension(2),intent(out):: xm
239 logical ,intent(out):: ff
240 integer(spi):: i
241 do i=1,2; call zttozm(a,xt(i),xm(i),ff); if(ff)return; enddo
242 end subroutine xttoxm
243 
253 subroutine xmtoxt(a,xm,xt,xtd,ff)! [xmtoxt]
254 implicit none
255 real(dp), intent(in ):: a
256 real(dp),dimension(2), intent(in ):: xm
257 real(dp),dimension(2), intent(out):: xt
258 real(dp),dimension(2,2),intent(out):: xtd
259 logical, intent(out):: ff
260 integer(spi):: i
261 xtd=u0; do i=1,2; call zmtozt(a,xm(i),xt(i),xtd(i,i),ff); if(ff)return; enddo
262 end subroutine xmtoxt
263 
275 subroutine xmtoxt1(a,xm,xt,xtd,xt1,xtd1,ff)! [xmtoxt]
276 implicit none
277 real(dp), intent(in ):: a
278 real(dp),dimension(2), intent(in ):: xm
279 real(dp),dimension(2), intent(out):: xt,xt1
280 real(dp),dimension(2,2), intent(out):: xtd,xtd1
281 logical, intent(out):: ff
282 integer(spi):: i
283 xtd=u0
284 xtd1=u0
285 do i=1,2
286  call zmtozt1(a,xm(i),xt(i),xtd(i,i),xt1(i),xtd1(i,i),ff)
287  if(ff)return
288 enddo
289 end subroutine xmtoxt1
290 
298 subroutine zttozm(a,zt,zm,ff)! [zttozm]
299 implicit none
300 real(dp),intent(in ):: a,zt
301 real(dp),intent(out):: zm
302 logical, intent(out):: ff
303 real(dp):: ra,razt
304 ff=f
305 if (a>u0)then; ra=sqrt( a); razt=ra*zt; zm=atan(razt)/ra
306 elseif(a<u0)then; ra=sqrt(-a); razt=ra*zt; ff=abs(razt)>=u1; if(ff)return
307  zm=atanh(razt)/ra
308 else ; zm=zt
309 endif
310 end subroutine zttozm
311 
322 subroutine zmtozt(a,zm,zt,ztd,ff)! [zmtozt]
323 implicit none
324 real(dp),intent(in ):: a,zm
325 real(dp),intent(out):: zt,ztd
326 logical, intent(out):: ff
327 real(dp):: ra
328 ff=f
329 if (a>u0)then; ra=sqrt( a); zt=tan(ra*zm)/ra; ff=abs(ra*zm)>=pih
330 elseif(a<u0)then; ra=sqrt(-a); zt=tanh(ra*zm)/ra
331 else ; zt=zm
332 endif
333 ztd=u1+a*zt*zt
334 end subroutine zmtozt
335 
347 subroutine zmtozt1(a,zm,zt,ztd,zt1,ztd1,ff)! [zmtozt]
348 use pietc, only: o3
349 use pfun, only: sinoxm,sinox,sinhoxm,sinhox
350 implicit none
351 real(dp),intent(in ):: a,zm
352 real(dp),intent(out):: zt,ztd,zt1,ztd1
353 logical, intent(out):: ff
354 real(dp):: ra,rad,razm
355 ff=f
356 if (a>u0)then;ra=sqrt( a);razm=ra*zm; zt=tan(razm)/ra; ff=abs(razm)>=pih
357 rad=o2/ra
358 zt1=(rad*zm/ra)*((-u2*sin(razm*o2)**2-sinoxm(razm))/cos(razm)+(tan(razm))**2)
359 elseif(a<u0)then;ra=sqrt(-a);razm=ra*zm; zt=tanh(razm)/ra
360 rad=-o2/ra
361 zt1=(rad*zm/ra)*((u2*sinh(razm*o2)**2-sinhoxm(razm))/cosh(razm)-(tanh(razm))**2)
362 else ;zt=zm; zt1=zm**3*o3
363 endif
364 ztd=u1+a*zt*zt
365 ztd1=zt*zt +u2*a*zt*zt1
366 end subroutine zmtozt1
367 
380 subroutine xmtoxc_vak(ak,xm,xc,xcd,ff)! [xmtoxc_ak]
381 implicit none
382 real(dp),dimension(2), intent(in ):: ak,xm
383 real(dp),dimension(3), intent(out):: xc
384 real(dp),dimension(3,2),intent(out):: xcd
385 logical, intent(out):: ff
386 call xmtoxc_ak(ak(1),ak(2),xm,xc,xcd,ff)
387 end subroutine xmtoxc_vak
388 
399 subroutine xmtoxc_vak1(ak,xm,xc,xcd,xc1,xcd1,ff)! [xmtoxc_ak]
400 implicit none
401 real(dp),dimension(2), intent(in ):: ak,xm
402 real(dp),dimension(3), intent(out):: xc
403 real(dp),dimension(3,2), intent(out):: xcd
404 real(dp),dimension(3,2), intent(out):: xc1
405 real(dp),dimension(3,2,2),intent(out):: xcd1
406 logical, intent(out):: ff
407 real(dp),dimension(3,2,2):: xcdd
408 real(dp),dimension(2,2,2):: xsd1,xsdd
409 real(dp),dimension(2,2) :: xtd,xsd,xs1,xtd1,xsdk
410 real(dp),dimension(2) :: xt,xt1,xs,xsk
411 integer(spi) :: i
412 call xmtoxt1(ak(1),xm,xt,xtd,xt1,xtd1,ff); if(ff)return
413 call xttoxs1(ak(2),xt,xs,xsd,xsdd,xsk,xsdk,ff); if(ff)return
414 xs1(:,2)=xsk; xs1(:,1)=matmul(xsd,xt1)
415 xsd1(:,:,1)=matmul(xsd,xtd1)+matmul(xsdd(:,:,1)*xt1(1)+xsdd(:,:,2)*xt1(2),xtd)
416 xsd1(:,:,2)=matmul(xsdk,xtd)
417 xsd=matmul(xsd,xtd)
418 call xstoxc(xs,xc,xcd,xcdd)
419 xc1=matmul(xcd,xs1)
420 do i=1,3; xcd1(i,:,:)=matmul(transpose(xsd),matmul(xcdd(i,:,:),xs1)); enddo
421 do i=1,2; xcd1(:,:,i)=xcd1(:,:,i)+matmul(xcd,xsd1(:,:,i)); enddo
422 xcd=matmul(xcd,xsd)
423 end subroutine xmtoxc_vak1
424 
438 subroutine xmtoxc_ak(a,k,xm,xc,xcd,ff)! [xmtoxc_ak]
439 implicit none
440 real(dp), intent(in ):: a,k
441 real(dp),dimension(2), intent(in ):: xm
442 real(dp),dimension(3), intent(out):: xc
443 real(dp),dimension(3,2),intent(out):: xcd
444 logical, intent(out):: ff
445 real(dp),dimension(2,2):: xtd,xsd
446 real(dp),dimension(2) :: xt,xs
447 call xmtoxt(a,xm,xt,xtd,ff); if(ff)return
448 call xttoxs(k,xt,xs,xsd,ff); if(ff)return
449 xsd=matmul(xsd,xtd)
450 call xstoxc(xs,xc,xcd)
451 xcd=matmul(xcd,xsd)
452 end subroutine xmtoxc_ak
453 
464 subroutine xctoxm_ak(a,k,xc,xm,ff)! [xctoxm_ak]
465 implicit none
466 real(dp), intent(in ):: a,k
467 real(dp),dimension(3),intent(in ):: xc
468 real(dp),dimension(2),intent(out):: xm
469 logical, intent(out):: ff
470 real(dp),dimension(2):: xs,xt
471 ff=f
472 call xctoxs(xc,xs)
473 call xstoxt(k,xs,xt,ff); if(ff)return
474 call xttoxm(a,xt,xm,ff)
475 end subroutine xctoxm_ak
476 
487 subroutine get_edges(arcx,arcy,edgex,edgey)! [get_edges]
488 implicit none
489 real(dp), intent(in ):: arcx,arcy
490 real(dp),dimension(3),intent(out):: edgex,edgey
491 real(dp):: cx,sx,cy,sy,rarcx,rarcy
492 rarcx=arcx*dtor; rarcy=arcy*dtor
493 cx=cos(rarcx); sx=sin(rarcx)
494 cy=cos(rarcy); sy=sin(rarcy)
495 edgex=(/sx,u0,cx/); edgey=(/u0,sy,cy/)
496 end subroutine get_edges
497 
510 subroutine get_qx(j0, v1,v2,v3,v4)! [get_qx]
511 use psym2, only: logsym2
512 implicit none
513 real(dp),dimension(3,2),intent(in ):: j0
514 real(dp), intent(out):: v1,v2,v3,v4
515 real(dp),dimension(2,2):: el,g
516 g=matmul(transpose(j0),j0)
517 call logsym2(g,el); el=el*o2
518 v1=el(1,1)**2+u2*el(1,2)**2+el(2,2)**2
519 v2=el(1,1)
520 v3=el(2,2)
521 v4=(el(1,1)+el(2,2))**2
522 end subroutine get_qx
523 
540 subroutine get_qxd(j0,j0d, v1,v2,v3,v4,v1d,v2d,v3d,v4d)! [get_qx]
541 use psym2, only: logsym2
542 implicit none
543 real(dp),dimension(3,2), intent(in ):: j0
544 real(dp),dimension(3,2,2),intent(in ):: j0d
545 real(dp), intent(out):: v1,v2,v3,v4
546 real(dp),dimension(2), intent(out):: v1d,v2d,v3d,v4d
547 real(dp),dimension(2,2,2,2):: deldg
548 real(dp),dimension(2,2,2) :: eld,gd
549 real(dp),dimension(2,2) :: el,g
550 integer(spi) :: i,j,k
551 g=matmul(transpose(j0),j0)
552 do i=1,2
553  gd(:,:,i)=matmul(transpose(j0d(:,:,i)),j0)+matmul(transpose(j0),j0d(:,:,i))
554 enddo
555 call logsym2(g,el,deldg); el=el*o2; deldg=deldg*o2
556 eld=u0
557 do i=1,2; do j=1,2; do k=1,2
558  eld(:,:,k)=eld(:,:,k)+deldg(:,:,i,j)*gd(i,j,k)
559 enddo ; enddo ; enddo
560 v1=el(1,1)**2+u2*el(1,2)**2+el(2,2)**2
561 v2=el(1,1)
562 v3=el(2,2)
563 v4=(el(1,1)+el(2,2))**2
564 v1d=u2*(el(1,1)*eld(1,1,:)+u2*el(1,2)*eld(1,2,:)+el(2,2)*eld(2,2,:))
565 v2d=eld(1,1,:)
566 v3d=eld(2,2,:)
567 v4d=u2*(el(1,1)+el(2,2))*(eld(1,1,:)+eld(2,2,:))
568 end subroutine get_qxd
569 
601 subroutine get_meanqd(ngh,lam,xg,wg,ak,ma, q,qdak,qdma, & ! [get_meanq]
602  ga,gadak,gadma, ff)
603 implicit none
604 integer(spi), intent(in ):: ngh
605 real(dp), intent(in ):: lam
606 real(dp),dimension(ngh),intent(in ):: xg,wg
607 real(dp),dimension(2) ,intent(in ):: ak,ma
608 real(dp), intent(out):: q
609 real(dp),dimension(2), intent(out):: qdak,qdma
610 real(dp),dimension(2), intent(out):: ga
611 real(dp),dimension(2,2),intent(out):: gadak,gadma
612 logical, intent(out):: ff
613 real(dp),dimension(3,2,2):: xcd1
614 real(dp),dimension(3,2) :: xcd,xc1
615 real(dp),dimension(3) :: xc
616 real(dp),dimension(2) :: xm, v1dxy,v2dxy,v3dxy,v4dxy, &
617  v1dl,v2dl,v3dl,v4dl, v1d,v2d,v3d,v4d
618 real(dp) :: wx,wy, &
619  v1xy,v2xy,v3xy,v4xy, v1l,v2l,v3l,v4l, v1,v2,v3,v4
620 integer(spi) :: i,ic,ix,iy
621 v1 =u0; v2 =u0; v3 =u0; v4 =u0
622 v1d=u0; v2d=u0; v3d=u0; v4d=u0
623 do iy=1,ngh
624  wy=wg(iy)
625  xm(2)=ma(2)*xg(iy)
626  v1l =u0; v2l =u0; v3l =u0; v4l =u0
627  v1dl=u0; v2dl=u0; v3dl=u0; v4dl=u0
628  do ix=1,ngh
629  wx=wg(ix)
630  xm(1)=ma(1)*xg(ix)
631  call xmtoxc_ak(ak,xm,xc,xcd,xc1,xcd1,ff); if(ff)return
632  call get_qx(xcd,xcd1, v1xy,v2xy,v3xy,v4xy, v1dxy,v2dxy,v3dxy,v4dxy)
633  v1l =v1l +wx*v1xy; v2l =v2l +wx*v2xy
634  v3l =v3l +wx*v3xy; v4l =v4l +wx*v4xy
635  v1dl=v1dl+wx*v1dxy;v2dl=v2dl+wx*v2dxy
636  v3dl=v3dl+wx*v3dxy;v4dl=v4dl+wx*v4dxy
637  enddo
638  v1 =v1 +wy*v1l; v2 =v2 +wy*v2l; v3 =v3 +wy*v3l; v4 =v4 +wy*v4l
639  v1d=v1d+wy*v1dl; v2d=v2d+wy*v2dl; v3d=v3d+wy*v3dl; v4d=v4d+wy*v4dl
640 enddo
641 call get_qofv(lam,v1,v2,v3,v4, q)! <- Q(lam) based on the v1,v2,v3,v4
642 call get_qofv(lam,v2,v3, v1d,v2d,v3d,v4d, qdak)! <- Derivative of Q wrt ak
643 ! Derivatives of ga wrt ak, and of q and ga wrt ma:
644 gadma=u0! <- needed because only diagonal elements are filled
645 do i=1,2
646  ic=3-i
647  xm=0; xm(i)=ma(i)
648  call xmtoxc_ak(ak,xm,xc,xcd,xc1,xcd1,ff); if(ff)return
649  ga(i)=atan2(xc(i),xc(3))*rtod
650  gadak(i,:)=(xc(3)*xc1(i,:)-xc(i)*xc1(3,:))*rtod
651  gadma(i,i)=(xc(3)*xcd(i,i)-xc(i)*xcd(3,i))*rtod
652 
653  v1l=u0; v2l=u0; v3l=u0; v4l=u0
654  do iy=1,ngh
655  wy=wg(iy)
656  xm(ic)=ma(ic)*xg(iy)
657  call xmtoxc_ak(ak,xm,xc,xcd,ff); if(ff)return
658  call get_qx(xcd, v1xy,v2xy,v3xy,v4xy)
659  v1l=v1l+wy*v1xy; v2l=v2l+wy*v2xy; v3l=v3l+wy*v3xy; v4l=v4l+wy*v4xy
660  enddo
661  v1d(i)=(v1l-v1)/ma(i); v2d(i)=(v2l-v2)/ma(i)
662  v3d(i)=(v3l-v3)/ma(i); v4d(i)=(v4l-v4)/ma(i)
663 enddo
664 call get_qofv(lam,v2,v3, v1d,v2d,v3d,v4d, qdma)
665 end subroutine get_meanqd
666 
680 subroutine get_meanqs(n,ngh,lam,xg,wg,aks,mas, qs,ff)! [get_meanq]
681 implicit none
682 integer(spi), intent(in ):: n,ngh
683 real(dp),dimension(ngh),intent(in ):: xg,wg
684 real(dp), intent(in ):: lam
685 real(dp),dimension(2,n),intent(in ):: aks,mas
686 real(dp),dimension(n), intent(out):: qs
687 logical, intent(out):: ff
688 real(dp),dimension(n) :: v1s,v2s,v3s,v4s
689 real(dp),dimension(n) :: v1sl,v2sl,v3sl,v4sl
690 real(dp),dimension(3,2):: xcd
691 real(dp),dimension(3) :: xc
692 real(dp),dimension(2) :: xm,xgs
693 real(dp) :: wx,wy, v1xy,v2xy,v3xy,v4xy
694 integer(spi) :: i,ix,iy
695 v1s=u0; v2s=u0; v3s=u0; v4s=u0
696 do iy=1,ngh
697  wy=wg(iy)
698  v1sl=u0; v2sl=u0; v3sl=u0; v4sl=u0
699  do ix=1,ngh
700  wx=wg(ix)
701  xgs=(/xg(ix),xg(iy)/)
702  do i=1,n
703  xm=mas(:,i)*xgs
704  call xmtoxc_ak(aks(:,i),xm,xc,xcd,ff); if(ff)return
705  call get_qx(xcd,v1xy,v2xy,v3xy,v4xy)
706  v1sl(i)=v1sl(i)+wx*v1xy; v2sl(i)=v2sl(i)+wx*v2xy
707  v3sl(i)=v3sl(i)+wx*v3xy; v4sl(i)=v4sl(i)+wx*v4xy
708  enddo
709  enddo
710  v1s=v1s+wy*v1sl; v2s=v2s+wy*v2sl; v3s=v3s+wy*v3sl; v4s=v4s+wy*v4sl
711 enddo
712 call get_qofv(n,lam,v1s,v2s,v3s,v4s, qs)
713 end subroutine get_meanqs
714 
727 subroutine get_qofv(lam,v1,v2,v3,v4, q)! [get_qofv]
728 implicit none
729 real(dp),intent(in ):: lam,v1,v2,v3,v4
730 real(dp),intent(out):: q
731 real(dp):: lamc
732 lamc=u1-lam
733 q=lamc*(v1-(v2**2+v3**2)) +lam*(v4 -(v2+v3)**2)
734 end subroutine get_qofv
735 
749 subroutine get_qofvd(lam, v2,v3, v1d,v2d,v3d,v4d, qd)! [get_qofv]
750 implicit none
751 real(dp), intent(in ):: lam,v2,v3
752 real(dp),dimension(2),intent(in ):: v1d,v2d,v3d,v4d
753 real(dp),dimension(2),intent(out):: qd
754 real(dp):: lamc
755 lamc=u1-lam
756 qd=lamc*(v1d-u2*(v2d*v2+v3d*v3))+lam*(v4d-u2*(v2d+v3d)*(v2+v3))
757 end subroutine get_qofvd
758 
769 subroutine get_qsofvs(n,lam,v1s,v2s,v3s,v4s, qs)! [get_qofv]
770 implicit none
771 integer(spi), intent(in ):: n
772 real(dp), intent(in ):: lam
773 real(dp),dimension(n),intent(in ):: v1s,v2s,v3s,v4s
774 real(dp),dimension(n),intent(out):: qs
775 real(dp):: lamc
776 lamc=u1-lam
777 qs=lamc*(v1s-(v2s**2+v3s**2)) +lam*(v4s -(v2s+v3s)**2)
778 end subroutine get_qsofvs
779 
789 subroutine guessak_map(asp,tmarcx,ak)! [guessak_map]
790 implicit none
791 real(dp), intent(in ):: asp,tmarcx
792 real(dp),dimension(2),intent(out):: ak
793 real(dp):: gmarcx
794 gmarcx=tmarcx*rtod
795 call guessak_geo(asp,gmarcx,ak)
796 end subroutine guessak_map
797 
807 subroutine guessak_geo(asp,arc,ak)! [guessak_geo]
808 implicit none
809 real(dp), intent(in ):: asp,arc
810 real(dp),dimension(2),intent(out):: ak
811 integer(spi),parameter :: narc=11,nasp=10! <- Table index bounds
812 real(dp), parameter :: eps=1.e-7_dp,darc=10._dp+eps,dasp=.1_dp+eps
813 real(dp),dimension(nasp,0:narc):: adarc,kdarc
814 real(dp) :: sasp,sarc,wx0,wx1,wa0,wa1
815 integer(spi) :: iasp0,iasp1,iarc0,iarc1
816 !------------------------
817 ! Tables of approximate A (adarc) and K (kdarc), valid for lam=0.8, for aspect ratio,
818 ! asp, at .1, .2, .3, .4, .5, .6, .7, .8, 1. and major semi-arc, arcx, at values,
819 ! 0 (nominally), 10., ..., 100. degrees (where the nominal "0" angle is actually
820 ! from a computation at 3 degrees, since zero would not make sense).
821 ! The 100 degree rows are repeated as the 110 degree entries deliberately to pad the
822 ! table into the partly forbidden parameter space to allow skinny domains of up to
823 ! 110 degree semi-length to be validly endowed with optimal A and K:
824 data adarc/ &
825 -.450_dp,-.328_dp,-.185_dp,-.059_dp,0.038_dp,0.107_dp,0.153_dp,0.180_dp,0.195_dp,0.199_dp,&
826 -.452_dp,-.327_dp,-.184_dp,-.058_dp,0.039_dp,0.108_dp,0.154_dp,0.182_dp,0.196_dp,0.200_dp,&
827 -.457_dp,-.327_dp,-.180_dp,-.054_dp,0.043_dp,0.112_dp,0.158_dp,0.186_dp,0.200_dp,0.205_dp,&
828 -.464_dp,-.323_dp,-.173_dp,-.047_dp,0.050_dp,0.118_dp,0.164_dp,0.192_dp,0.208_dp,0.213_dp,&
829 -.465_dp,-.313_dp,-.160_dp,-.035_dp,0.060_dp,0.127_dp,0.173_dp,0.202_dp,0.217_dp,0.224_dp,&
830 -.448_dp,-.288_dp,-.138_dp,-.017_dp,0.074_dp,0.140_dp,0.184_dp,0.213_dp,0.230_dp,0.237_dp,&
831 -.395_dp,-.244_dp,-.104_dp,0.008_dp,0.093_dp,0.156_dp,0.199_dp,0.227_dp,0.244_dp,0.252_dp,&
832 -.301_dp,-.177_dp,-.057_dp,0.042_dp,0.119_dp,0.175_dp,0.215_dp,0.242_dp,0.259_dp,0.269_dp,&
833 -.185_dp,-.094_dp,0.001_dp,0.084_dp,0.150_dp,0.199_dp,0.235_dp,0.260_dp,0.277_dp,0.287_dp,&
834 -.069_dp,-.006_dp,0.066_dp,0.132_dp,0.186_dp,0.227_dp,0.257_dp,0.280_dp,0.296_dp,0.308_dp,&
835 0.038_dp,0.081_dp,0.134_dp,0.185_dp,0.226_dp,0.258_dp,0.283_dp,0.303_dp,0.319_dp,0.333_dp,&
836 0.038_dp,0.081_dp,0.134_dp,0.185_dp,0.226_dp,0.258_dp,0.283_dp,0.303_dp,0.319_dp,0.333_dp/
837 
838 data kdarc/ &
839 -.947_dp,-.818_dp,-.668_dp,-.535_dp,-.433_dp,-.361_dp,-.313_dp,-.284_dp,-.269_dp,-.264_dp,&
840 -.946_dp,-.816_dp,-.665_dp,-.533_dp,-.431_dp,-.359_dp,-.311_dp,-.282_dp,-.267_dp,-.262_dp,&
841 -.942_dp,-.806_dp,-.655_dp,-.524_dp,-.424_dp,-.353_dp,-.305_dp,-.276_dp,-.261_dp,-.255_dp,&
842 -.932_dp,-.789_dp,-.637_dp,-.509_dp,-.412_dp,-.343_dp,-.296_dp,-.266_dp,-.250_dp,-.244_dp,&
843 -.909_dp,-.759_dp,-.609_dp,-.486_dp,-.394_dp,-.328_dp,-.283_dp,-.254_dp,-.237_dp,-.230_dp,&
844 -.863_dp,-.711_dp,-.569_dp,-.456_dp,-.372_dp,-.310_dp,-.266_dp,-.238_dp,-.220_dp,-.212_dp,&
845 -.779_dp,-.642_dp,-.518_dp,-.419_dp,-.343_dp,-.287_dp,-.247_dp,-.220_dp,-.202_dp,-.192_dp,&
846 -.661_dp,-.556_dp,-.456_dp,-.374_dp,-.310_dp,-.262_dp,-.226_dp,-.200_dp,-.182_dp,-.171_dp,&
847 -.533_dp,-.462_dp,-.388_dp,-.325_dp,-.274_dp,-.234_dp,-.203_dp,-.179_dp,-.161_dp,-.150_dp,&
848 -.418_dp,-.373_dp,-.322_dp,-.275_dp,-.236_dp,-.204_dp,-.178_dp,-.156_dp,-.139_dp,-.127_dp,&
849 -.324_dp,-.296_dp,-.261_dp,-.229_dp,-.200_dp,-.174_dp,-.152_dp,-.133_dp,-.117_dp,-.104_dp,&
850 -.324_dp,-.296_dp,-.261_dp,-.229_dp,-.200_dp,-.174_dp,-.152_dp,-.133_dp,-.117_dp,-.104_dp/
851 !=============================================================================
852 sasp=asp/dasp
853 iasp0=floor(sasp); wa1=sasp-iasp0
854 iasp1=iasp0+1; wa0=iasp1-sasp
855 sarc=arc/darc
856 iarc0=floor(sarc); wx1=sarc-iarc0
857 iarc1=iarc0+1; wx0=iarc1-sarc
858 if(iasp0<1 .or. iasp1>nasp)stop 'Guessak_geo; Aspect ratio out of range'
859 if(iarc0<0 .or. iarc1>narc)stop 'Guessak_geo; Major semi-arc is out of range'
860 
861 ! Bilinearly interpolate A and K from crude table into a 2-vector:
862 ak=(/wx0*(wa0*adarc(iasp0,iarc0)+wa1*adarc(iasp1,iarc0))+ &
863  wx1*(wa0*adarc(iasp0,iarc1)+wa1*adarc(iasp1,iarc1)), &
864  wx0*(wa0*kdarc(iasp0,iarc0)+wa1*kdarc(iasp1,iarc0))+ &
865  wx1*(wa0*kdarc(iasp0,iarc1)+wa1*kdarc(iasp1,iarc1))/)
866 end subroutine guessak_geo
867 
901 subroutine bestesg_geo(lam,garcx,garcy, a,k,marcx,marcy,q,ff)! [bestesg_geo]
902 use pietc, only: u5,o5,s18,s36,s54,s72,ms18,ms36,ms54,ms72
903 use pmat, only: inv
904 use psym2, only: chol2
905 implicit none
906 real(dp),intent(in ):: lam,garcx,garcy
907 real(dp),intent(out):: a,k,marcx,marcy,q
908 logical ,intent(out):: ff
909 integer(spi),parameter :: nit=200,mit=20,ngh=25
910 real(dp) ,parameter :: u2o5=u2*o5,&
911  f18=u2o5*s18,f36=u2o5*s36,f54=u2o5*s54,&
912  f72=u2o5*s72,mf18=-f18,mf36=-f36,mf54=-f54,&
913  mf72=-f72,& !<- (Fourier transform coefficients)
914  r=0.001_dp,rr=r*r,dang=pi2*o5,crit=1.e-14_dp
915 real(dp),dimension(ngh) :: wg,xg
916 real(dp),dimension(0:4,0:4):: em5 ! <- Fourier matrix for 5 points
917 real(dp),dimension(0:4) :: qs
918 real(dp),dimension(2,0:4) :: aks,mas
919 real(dp),dimension(2,2) :: basis0,basis,hess,el,gadak,gadma,madga,madak
920 real(dp),dimension(2) :: ak,dak,dma,vec2,grad,qdak,qdma,ga,ma,gat
921 real(dp) :: s,tgarcx,tgarcy,asp,ang
922 integer(spi) :: i,it
923 logical :: flip
924 data em5/o5,u2o5, u0,u2o5, u0,& ! <-The Fourier matrix for 5 points. Applied
925  o5, f18, f72,mf54, f36,& ! to the five 72-degree spaced values in a
926  o5,mf54, f36, f18,mf72,& ! column-vector, the product vector has the
927  o5,mf54,mf36, f18, f72,& ! components, wave-0, cos and sin wave-1,
928  o5, f18,mf72,mf54,mf36/ ! cos and sin wave-2.
929 ! First guess upper-triangular basis regularizing the samples used to
930 ! estimate the Hessian of q by finite differencing:
931 data basis0/0.1_dp,u0, 0.3_dp,0.3_dp/
932 ff=lam<u0 .or. lam>=u1
933 if(ff)then; print'("In bestesg_geo; lam out of range")';return; endif
934 ff= garcx<=u0 .or. garcy<=u0
935 if(ff)then
936  print'("In bestesg_geo; a nonpositive domain parameter, garcx or garcy")'
937  return
938 endif
939 call gaulegh(ngh,xg,wg)! <- Prepare Gauss-Legendre nodes xg and weights wg
940 flip=garcy>garcx
941 if(flip)then; tgarcx=garcy; tgarcy=garcx! <- Switch
942 else ; tgarcx=garcx; tgarcy=garcy! <- Don't switch
943 endif
944 ga=(/tgarcx,tgarcy/)
945 asp=tgarcy/tgarcx
946 basis=basis0
947 
948 call guessak_geo(asp,tgarcx,ak)
949 ma=ga*dtor*0.9_dp ! Shrink first estimate, to start always within bounds
950 
951 ! Perform a Newton iteration (except with imperfect Hessian) to find the
952 ! parameter vector, ak, at which the derivative of Q at constant geographical
953 ! semi-axes, ga, as given, goes to zero. The direct evaluation of the
954 ! Q-derivative at constant ma (which is what is actually computed in
955 ! get_meanq) therefore needs modification to obtain Q-derivarive at constant
956 ! ga:
957 ! dQ/d(ak)|_ga = dQ/d(ak)|_ma - dQ/d(ma)|_ak*d(ma)/d(ga)|_ak*d(ga)/d(ak)|_ma
958 !
959 ! Since the Hessian evaluation is only carried out at constant map-space
960 ! semi-axes, ma, it is not ideal for this problem; consequently, the allowance
961 ! of newton iterations, nit, is much more liberal than we allow for the
962 ! companion routine, bestesg_map, where the constant ma condition WAS
963 ! appropriate.
964 do it=1,nit
965  call get_meanq(ngh,lam,xg,wg,ak,ma,q,qdak,qdma,gat,gadak,gadma,ff)
966  if(ff)return
967  madga=gadma; call inv(madga)! <- d(ma)/d(ga)|_ak ("at constant ak")
968  madak=-matmul(madga,gadak)
969  qdak=qdak+matmul(qdma,madak)! dQ/d(ak)|_ga
970  if(it<=mit)then ! <- Only recompute aks if the basis is new
971 ! Place five additional sample points around the stencil-ellipse:
972  do i=0,4
973  ang=i*dang ! steps of 72 degrees
974  vec2=(/cos(ang),sin(ang)/)*r ! points on a circle of radius r ...
975  dak=matmul(basis,vec2)
976  dma=matmul(madak,dak)
977  aks(:,i)=ak+dak ! ... become points on an ellipse.
978  mas(:,i)=ma+dma
979  enddo
980  call get_meanq(5,ngh,lam,xg,wg,aks,mas, qs,ff)
981  endif
982  grad=matmul(qdak,basis)/q ! <- New grad estimate, accurate to near roundoff
983  if(it<mit)then ! Assume the hessian will have significantly changed
984 ! Recover an approximate Hessian estimate, normalized by q, from all
985 ! 6 samples, q, qs. These are wrt the basis, not (a,k) directly. Note that
986 ! this finite-difference Hessian is wrt q at constant ak, which is roughly
987 ! approximating the desired Hessian wrt q at constant ga; the disparity
988 ! is the reason that we allow more iterations in this routine than in
989 ! companion routine bestesg_map, since we cannot achieve superlinear
990 ! convergence in the present case.
991 ! Make qs the 5-pt discrete Fourier coefficients of the ellipse pts:
992  qs=matmul(em5,qs)/q
993 ! grad=qs(1:2)/r ! Old finite difference estimate is inferior to new grad
994  qs(0)=qs(0)-u1! <- cos(2*ang) coefficient relative to the central value.
995  hess(1,1)=qs(0)+qs(3)! <- combine cos(0*ang) and cos(2*ang) coefficients
996  hess(1,2)=qs(4) ! <- sin(2*ang) coefficient
997  hess(2,1)=qs(4) !
998  hess(2,2)=qs(0)-qs(3)! <- combine cos(0*ang) and cos(2*ang) coefficients
999  hess=hess*u2/rr ! <- rr is r**2
1000 
1001 ! Perform a Cholesky decomposition of the hessian:
1002  call chol2(hess,el,ff)
1003  if(ff)then
1004  print'(" In bestESG_geo, Hessian is not positive; cholesky fails")'
1005  return
1006  endif
1007 ! Invert the cholesky factor in place:
1008  el(1,1)=u1/el(1,1); el(2,2)=u1/el(2,2); el(2,1)=-el(2,1)*el(1,1)*el(2,2)
1009  endif
1010 ! Estimate a Newton step towards the minimum of function Q(a,k):
1011  vec2=-matmul(transpose(el),matmul(el,grad))
1012  dak=matmul(basis,vec2)
1013  gat=gat+matmul(gadak,dak)! <- increment ga
1014  ak=ak+dak! <- increment the parameter vector estimate
1015  dma=-matmul(madga,gat-ga)
1016  ma=ma+dma
1017 
1018 ! Use the inverse cholesky factor to re-condition the basis. This is to make
1019 ! the next stencil-ellipse more closely share the shape of the elliptical
1020 ! contours of Q near its minumum -- essentially a preconditioning of the
1021 ! numerical optimization:
1022  if(it<mit)basis=matmul(basis,transpose(el))
1023 
1024  s=sqrt(dot_product(dak,dak))
1025  if(s<crit)exit ! <-Sufficient convergence of the Newton iteration
1026 enddo
1027 if(it>nit)print'("WARNING; Relatively inferior convergence in bestesg_geo")'
1028 a=ak(1)
1029 k=ak(2)
1030 if(flip)then; marcx=ma(2); marcy=ma(1)! Remember to switch back
1031 else ; marcx=ma(1); marcy=ma(2)! don't switch
1032 endif
1033 end subroutine bestesg_geo
1034 
1071 subroutine bestesg_map(lam,marcx,marcy, a,k,garcx,garcy,q,ff) ![bestesg_map]
1072 use pietc, only: u5,o5,s18,s36,s54,s72,ms18,ms36,ms54,ms72
1073 use psym2, only: chol2
1074 implicit none
1075 real(dp),intent(in ):: lam,marcx,marcy
1076 real(dp),intent(out):: a,k,garcx,garcy,q
1077 logical ,intent(out):: ff
1078 integer(spi),parameter :: nit=25,mit=7,ngh=25
1079 real(dp),parameter :: u2o5=u2*o5, &
1080  f18=u2o5*s18,f36=u2o5*s36,f54=u2o5*s54, &
1081  f72=u2o5*s72,mf18=-f18,mf36=-f36,mf54=-f54,&
1082  mf72=-f72,& !<- (Fourier)
1083  r=0.001_dp,rr=r*r,dang=pi2*o5,crit=1.e-12_dp
1084 real(dp),dimension(ngh) :: wg,xg
1085 real(dp),dimension(0:4,0:4):: em5 ! <- Fourier matrix for 5 points
1086 real(dp),dimension(0:4) :: qs ! <-Sampled q, its Fourier coefficients
1087 real(dp),dimension(2,0:4) :: aks,mas! <- tiny elliptical array of ak
1088 real(dp),dimension(2,2) :: basis0,basis,hess,el,gadak,gadma
1089 real(dp),dimension(2) :: ak,dak,vec2,grad,qdak,qdma,ga,ma
1090 real(dp) :: s,tmarcx,tmarcy,asp,ang
1091 integer(spi) :: i,it
1092 logical :: flip
1093 data em5/o5,u2o5, u0,u2o5, u0,& ! <-The Fourier matrix for 5 points. Applied
1094  o5, f18, f72,mf54, f36,& ! to the five 72-degree spaced values in a
1095  o5,mf54, f36, f18,mf72,& ! column-vector, the product vector has the
1096  o5,mf54,mf36, f18, f72,& ! components, wave-0, cos and sin wave-1,
1097  o5, f18,mf72,mf54,mf36/ ! cos and sin wave-2.
1098 ! First guess upper-triangular basis regularizing the samples used to
1099 ! estimate the Hessian of q by finite differencing:
1100 data basis0/0.1_dp,u0, 0.3_dp,0.3_dp/
1101 ff=lam<u0 .or. lam>=u1
1102 if(ff)then; print'("In bestesg_map; lam out of range")';return; endif
1103 ff= marcx<=u0 .or. marcy<=u0
1104 if(ff)then
1105  print'("In bestesg_map; a nonpositive domain parameter, marcx or marcy")'
1106  return
1107 endif
1108 call gaulegh(ngh,xg,wg)
1109 flip=marcy>marcx
1110 if(flip)then; tmarcx=marcy; tmarcy=marcx! <- Switch
1111 else ; tmarcx=marcx; tmarcy=marcy! <- Don't switch
1112 endif
1113 ma=(/tmarcx,tmarcy/); do i=0,4; mas(:,i)=ma; enddo
1114 asp=tmarcy/tmarcx
1115 basis=basis0
1116 
1117 call guessak_map(asp,tmarcx,ak)
1118 
1119 do it=1,nit
1120  call get_meanq(ngh,lam,xg,wg,ak,ma,q,qdak,qdma,ga,gadak,gadma,ff)
1121  if(ff)return
1122  if(it<=mit)then
1123 ! Place five additional sample points around the stencil-ellipse:
1124  do i=0,4
1125  ang=i*dang ! steps of 72 degrees
1126  vec2=(/cos(ang),sin(ang)/)*r ! points on a circle of radius r ...
1127  aks(:,i)=ak+matmul(basis,vec2) ! ... become points on an ellipse.
1128  enddo
1129  call get_meanq(5,ngh,lam,xg,wg,aks,mas, qs,ff)
1130  endif
1131  grad=matmul(qdak,basis)/q ! <- New grad estimate, accurate to near roundoff
1132  if(it<mit)then
1133 ! Recover Hessian estimate, normalized by q, from all
1134 ! 6 samples, q, qs. These are wrt the basis, not (a,k) directly.
1135 ! The Hessian estimate uses a careful finite method, which is accurate
1136 ! enough. The gradient is NOT estimated by finite differences because we
1137 ! need the gradient to be accurate to near roundoff levels in order that
1138 ! the converged Newton iteration is a precise solution.
1139 ! Make qs the 5-pt discrete Fourier coefficients of the ellipse pts:
1140  qs=matmul(em5,qs)/q
1141 ! grad=qs(1:2)/r ! Old finite difference estimate is inferior to new grad
1142  qs(0)=qs(0)-u1
1143  hess(1,1)=qs(0)+qs(3)! <- combine cos(0*ang) and cos(2*ang) coefficients
1144  hess(1,2)=qs(4) ! <- sin(2*ang) coefficient
1145  hess(2,1)=qs(4) !
1146  hess(2,2)=qs(0)-qs(3)! <- combine cos(0*ang) and cos(2*ang) coefficients
1147  hess=hess*u2/rr ! <- rr is r**2
1148 
1149 ! Perform a Cholesky decomposition of the hessian:
1150  call chol2(hess,el,ff)
1151  if(ff)then
1152  print'(" In bestESG_map, hessian is not positive; cholesky fails")'
1153  return
1154  endif
1155 ! Invert the cholesky factor in place:
1156  el(1,1)=u1/el(1,1); el(2,2)=u1/el(2,2); el(2,1)=-el(2,1)*el(1,1)*el(2,2)
1157  endif
1158 ! Estimate a Newton step towards the minimum of function Q(a,k):
1159  vec2=-matmul(transpose(el),matmul(el,grad))
1160  dak=matmul(basis,vec2)
1161  ga=ga+matmul(gadak,dak)! <- increment ga
1162  ak=ak+dak! <- increment the parameter vector estimate
1163 
1164 ! Use the inverse cholesky factor to re-condition the basis. This is to make
1165 ! the next stencil-ellipse more closely share the shape of the elliptical
1166 ! contours of Q near its minumum -- essentially a preconditioning of the
1167 ! numerical optimization:
1168  if(it<mit)basis=matmul(basis,transpose(el))
1169 
1170  s=sqrt(dot_product(dak,dak))
1171  if(s<crit)exit ! <-Sufficient convergence of the Newton iteration
1172 enddo
1173 if(it>nit)print'("WARNING; Relatively inferior convergence in bestesg_map")'
1174 a=ak(1)
1175 k=ak(2)
1176 if(flip)then; garcx=ga(2); garcy=ga(1)! Remember to switch back
1177 else ; garcx=ga(1); garcy=ga(2)! don't switch
1178 endif
1179 end subroutine bestesg_map
1180 
1218 subroutine hgrid_ak_rr(lx,ly,nx,ny,A,K,plat,plon,pazi, & ! [hgrid_ak_rr]
1219  delx,dely, glat,glon,garea, ff)
1220 use pmat4, only: sarea
1221 use pmat5, only: ctogr
1222 implicit none
1223 integer(spi), intent(in ):: lx,ly,nx,ny
1224 real(dp), intent(in ):: a,k,plat,plon,pazi, &
1225  delx,dely
1226 real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: glat,glon
1227 real(dp),dimension(lx:lx+nx-1,ly:ly+ny-1),intent(out):: garea
1228 logical, intent(out):: ff
1229 real(dp),dimension(3,3):: prot,azirot
1230 real(dp),dimension(3,2):: xcd
1231 real(dp),dimension(3) :: xc
1232 real(dp),dimension(2) :: xm
1233 real(dp) :: clat,slat,clon,slon,cazi,sazi,&
1234  rlat,drlata,drlatb,drlatc, &
1235  rlon,drlona,drlonb,drlonc
1236 integer(spi) :: ix,iy,mx,my
1237 clat=cos(plat); slat=sin(plat)
1238 clon=cos(plon); slon=sin(plon)
1239 cazi=cos(pazi); sazi=sin(pazi)
1240 
1241 azirot(:,1)=(/ cazi, sazi, u0/)
1242 azirot(:,2)=(/-sazi, cazi, u0/)
1243 azirot(:,3)=(/ u0, u0, u1/)
1244 
1245 prot(:,1)=(/ -slon, clon, u0/)
1246 prot(:,2)=(/-slat*clon, -slat*slon, clat/)
1247 prot(:,3)=(/ clat*clon, clat*slon, slat/)
1248 prot=matmul(prot,azirot)
1249 mx=lx+nx ! Index of the 'right' edge of the rectangular grid
1250 my=ly+ny ! Index of the 'top' edge of the rectangular grid
1251 do iy=ly,my
1252  xm(2)=iy*dely
1253  do ix=lx,mx
1254  xm(1)=ix*delx
1255  call xmtoxc_ak(a,k,xm,xc,xcd,ff)
1256  if(ff)return
1257  xcd=matmul(prot,xcd)
1258  xc =matmul(prot,xc )
1259  call ctogr(xc,glat(ix,iy),glon(ix,iy))
1260  enddo
1261 enddo
1262 
1263 ! Compute the areas of the quadrilateral grid cells:
1264 do iy=ly,my-1
1265  do ix=lx,mx-1
1266  rlat =glat(ix ,iy )
1267  drlata=glat(ix+1,iy )-rlat
1268  drlatb=glat(ix+1,iy+1)-rlat
1269  drlatc=glat(ix ,iy+1)-rlat
1270  rlon =glon(ix ,iy )
1271  drlona=glon(ix+1,iy )-rlon
1272  drlonb=glon(ix+1,iy+1)-rlon
1273  drlonc=glon(ix ,iy+1)-rlon
1274 ! If 'I' is the grid point (ix,iy), 'A' is (ix+1,iy); 'B' is (ix+1,iy+1)
1275 ! and 'C' is (ix,iy+1), then the area of the grid cell IABC is the sum of
1276 ! the areas of the traingles, IAB and IBC (the latter being the negative
1277 ! of the signed area of triangle, ICB):
1278  garea(ix,iy)=sarea(rlat, drlata,drlona, drlatb,drlonb) &
1279  -sarea(rlat, drlatc,drlonc, drlatb,drlonb)
1280  enddo
1281 enddo
1282 end subroutine hgrid_ak_rr
1283 
1335 subroutine hgrid_ak_rr_c(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak_rr]
1336  delx,dely, glat,glon,garea,dx,dy,angle_dx,angle_dy, ff)
1337 use pmat4, only: cross_product,triple_product
1338 use pmat5, only: ctogr
1339 implicit none
1340 integer(spi), intent(in ):: lx,ly,nx,ny
1341 real(dp), intent(in ):: a,k,plat,plon,pazi, &
1342  delx,dely
1343 real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: glat,glon
1344 real(dp),dimension(lx:lx+nx-1,ly:ly+ny-1),intent(out):: garea
1345 real(dp),dimension(lx:lx+nx-1,ly:ly+ny ),intent(out):: dx
1346 real(dp),dimension(lx:lx+nx ,ly:ly+ny-1),intent(out):: dy
1347 real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: angle_dx,angle_dy
1348 logical, intent(out):: ff
1349 real(dp),dimension(lx-1:lx+nx+1,ly-1:ly+ny+1):: gat ! Temporary area array
1350 real(dp),dimension(lx-1:lx+nx+1,ly :ly+ny ):: dxt ! Temporary dx array
1351 real(dp),dimension(lx :lx+nx ,ly-1:ly+ny+1):: dyt ! Temporary dy array
1352 real(dp),dimension(3,3):: prot,azirot
1353 real(dp),dimension(3,2):: xcd,eano
1354 real(dp),dimension(2,2):: xcd2
1355 real(dp),dimension(3) :: xc,east,north
1356 real(dp),dimension(2) :: xm
1357 real(dp) :: clat,slat,clon,slon,cazi,sazi,delxy
1358 integer(spi) :: ix,iy,mx,my,lxm,lym,mxp,myp
1359 delxy=delx*dely
1360 clat=cos(plat); slat=sin(plat)
1361 clon=cos(plon); slon=sin(plon)
1362 cazi=cos(pazi); sazi=sin(pazi)
1363 
1364 azirot(:,1)=(/ cazi, sazi, u0/)
1365 azirot(:,2)=(/-sazi, cazi, u0/)
1366 azirot(:,3)=(/ u0, u0, u1/)
1367 
1368 prot(:,1)=(/ -slon, clon, u0/)
1369 prot(:,2)=(/-slat*clon, -slat*slon, clat/)
1370 prot(:,3)=(/ clat*clon, clat*slon, slat/)
1371 prot=matmul(prot,azirot)
1372 
1373 mx=lx+nx ! Index of the 'right' edge of the rectangular grid
1374 my=ly+ny ! Index of the 'top' edge of the rectangular grid
1375 lxm=lx-1; mxp=mx+1 ! Indices of extra left and right edges
1376 lym=ly-1; myp=my+1 ! Indices of extra bottom and top edges
1377 
1378 !-- main body of horizontal grid:
1379 do iy=ly,my
1380  xm(2)=iy*dely
1381  do ix=lx,mx
1382  xm(1)=ix*delx
1383  call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
1384  xcd=matmul(prot,xcd)
1385  xc =matmul(prot,xc )
1386  call ctogr(xc,glat(ix,iy),glon(ix,iy))
1387  east=(/-xc(2),xc(1),u0/); east=east/sqrt(dot_product(east,east))
1388  north=cross_product(xc,east)
1389  eano(:,1)=east; eano(:,2)=north
1390  xcd2=matmul(transpose(eano),xcd)
1391  angle_dx(ix,iy)=atan2( xcd2(2,1),xcd2(1,1))
1392  angle_dy(ix,iy)=atan2(-xcd2(1,2),xcd2(2,2))
1393  dxt(ix,iy)=sqrt(dot_product(xcd2(:,1),xcd2(:,1)))*delx
1394  dyt(ix,iy)=sqrt(dot_product(xcd2(:,2),xcd2(:,2)))*dely
1395  gat(ix,iy)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
1396  enddo
1397 enddo
1398 
1399 !-- extra left edge, gat, dxt only:
1400 xm(1)=lxm*delx
1401 do iy=ly,my
1402  xm(2)=iy*dely
1403  call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
1404  xcd=matmul(prot,xcd)
1405  xc =matmul(prot,xc )
1406  east=(/-xc(2),xc(1),u0/); east=east/sqrt(dot_product(east,east))
1407  north=cross_product(xc,east)
1408  eano(:,1)=east; eano(:,2)=north
1409  xcd2=matmul(transpose(eano),xcd)
1410  dxt(lxm,iy)=sqrt(dot_product(xcd2(:,1),xcd2(:,1)))*delx
1411  gat(lxm,iy)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
1412 enddo
1413 
1414 !-- extra right edge, gat, dxt only:
1415 xm(1)=mxp*delx
1416 do iy=ly,my
1417  xm(2)=iy*dely
1418  call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
1419  xcd=matmul(prot,xcd)
1420  xc =matmul(prot,xc )
1421  east=(/-xc(2),xc(1),u0/); east=east/sqrt(dot_product(east,east))
1422  north=cross_product(xc,east)
1423  eano(:,1)=east; eano(:,2)=north
1424  xcd2=matmul(transpose(eano),xcd)
1425  dxt(mxp,iy)=sqrt(dot_product(xcd2(:,1),xcd2(:,1)))*delx
1426  gat(mxp,iy)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
1427 enddo
1428 
1429 !-- extra bottom edge, gat, dyt only:
1430 xm(2)=lym*dely
1431 do ix=lx,mx
1432  xm(1)=ix*delx
1433  call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
1434  xcd=matmul(prot,xcd)
1435  xc =matmul(prot,xc )
1436  east=(/-xc(2),xc(1),u0/); east=east/sqrt(dot_product(east,east))
1437  north=cross_product(xc,east)
1438  eano(:,1)=east; eano(:,2)=north
1439  xcd2=matmul(transpose(eano),xcd)
1440  dyt(ix,lym)=sqrt(dot_product(xcd2(:,2),xcd2(:,2)))*dely
1441  gat(ix,lym)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
1442 enddo
1443 
1444 !-- extra top edge, gat, dyt only:
1445 xm(2)=myp*dely
1446 do ix=lx,mx
1447  xm(1)=ix*delx
1448  call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
1449  xcd=matmul(prot,xcd)
1450  xc =matmul(prot,xc )
1451  east=(/-xc(2),xc(1),u0/); east=east/sqrt(dot_product(east,east))
1452  north=cross_product(xc,east)
1453  eano(:,1)=east; eano(:,2)=north
1454  xcd2=matmul(transpose(eano),xcd)
1455  dyt(ix,myp)=sqrt(dot_product(xcd2(:,2),xcd2(:,2)))*dely
1456  gat(ix,myp)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
1457 enddo
1458 
1459 ! Extra four corners, gat only:
1460 xm(2)=lym*dely
1461 !-- extra bottom left corner:
1462 xm(1)=lxm*delx
1463 call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
1464 xcd=matmul(prot,xcd)
1465 xc =matmul(prot,xc )
1466 gat(lxm,lym)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
1467 
1468 !-- extra bottom right corner:
1469 xm(1)=mxp*delx
1470 call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
1471 xcd=matmul(prot,xcd)
1472 xc =matmul(prot,xc )
1473 gat(mxp,lym)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
1474 
1475 xm(2)=myp*dely
1476 !-- extra top left corner:
1477 xm(1)=lxm*delx
1478 call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
1479 xcd=matmul(prot,xcd)
1480 xc =matmul(prot,xc )
1481 gat(lxm,myp)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
1482 
1483 !-- extra top right corner:
1484 xm(1)=mxp*delx
1485 call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
1486 xcd=matmul(prot,xcd)
1487 xc =matmul(prot,xc )
1488 gat(mxp,myp)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
1489 
1490 !-- 4th-order averaging over each central interval using 4-pt. stencils:
1491 dx =(13_dp*(dxt(lx :mx-1,:)+dxt(lx+1:mx ,:)) &
1492  -(dxt(lxm:mx-2,:)+dxt(lx+2:mxp,:)))/24_dp
1493 dy =(13_dp*(dyt(:,ly :my-1)+dyt(:,ly+1:my )) &
1494  -(dyt(:,lym:my-2)+dyt(:,ly+2:myp)))/24_dp
1495 gat(lx:mx-1,:)=(13_dp*(gat(lx :mx-1,:)+gat(lx+1:mx ,:)) &
1496  -(gat(lxm:mx-2,:)+gat(lx+2:mxp,:)))/24_dp
1497 garea =(13_dp*(gat(lx:mx-1,ly :my-1)+gat(lx:mx-1,ly+1:my )) &
1498  -(gat(lx:mx-1,lym:my-2)+gat(lx:mx-1,ly+2:myp)))/24_dp
1499 end subroutine hgrid_ak_rr_c
1500 
1538 subroutine hgrid_ak_rc(lx,ly,nx,ny,A,K,plat,plon,pazi, & ! [hgrid_ak_rc]
1539  delx,dely, xc,xcd,garea, ff)
1540 use pmat4, only: sarea
1541 use pmat5, only: ctogr
1542 implicit none
1543 integer(spi),intent(in ):: lx,ly,nx,ny
1544 real(dp), intent(in ):: a,k,plat,plon,pazi,delx,dely
1545 real(dp),dimension(3, lx:lx+nx ,ly:ly+ny ),intent(out):: xc
1546 real(dp),dimension(3,2,lx:lx+nx ,ly:ly+ny ),intent(out):: xcd
1547 real(dp),dimension( lx:lx+nx-1,ly:ly+ny-1),intent(out):: garea
1548 logical, intent(out):: ff
1549 real(dp),dimension(3,3):: prot,azirot
1550 real(dp),dimension(2) :: xm
1551 real(dp) :: clat,slat,clon,slon,cazi,sazi, &
1552  rlat,rlata,rlatb,rlatc,drlata,drlatb,drlatc, &
1553  rlon,rlona,rlonb,rlonc,drlona,drlonb,drlonc
1554 integer(spi) :: ix,iy,mx,my
1555 clat=cos(plat); slat=sin(plat)
1556 clon=cos(plon); slon=sin(plon)
1557 cazi=cos(pazi); sazi=sin(pazi)
1558 
1559 azirot(:,1)=(/ cazi, sazi, u0/)
1560 azirot(:,2)=(/-sazi, cazi, u0/)
1561 azirot(:,3)=(/ u0, u0, u1/)
1562 
1563 prot(:,1)=(/ -slon, clon, u0/)
1564 prot(:,2)=(/-slat*clon, -slat*slon, clat/)
1565 prot(:,3)=(/ clat*clon, clat*slon, slat/)
1566 prot=matmul(prot,azirot)
1567 mx=lx+nx ! Index of the 'right' edge of the rectangular grid
1568 my=ly+ny ! Index of the 'top' edge of the rectangular grid
1569 do iy=ly,my
1570  xm(2)=iy*dely
1571  do ix=lx,mx
1572  xm(1)=ix*delx
1573  call xmtoxc_ak(a,k,xm,xc(:,ix,iy),xcd(:,:,ix,iy),ff)
1574  if(ff)return
1575  xcd(:,:,ix,iy)=matmul(prot,xcd(:,:,ix,iy))
1576  xc(:, ix,iy)=matmul(prot,xc(:, ix,iy))
1577  enddo
1578 enddo
1579 
1580 ! Compute the areas of the quadrilateral grid cells:
1581 do iy=ly,my-1
1582  do ix=lx,mx-1
1583  call ctogr(xc(:,ix ,iy ),rlat ,rlon )
1584  call ctogr(xc(:,ix+1,iy ),rlata,rlona)
1585  call ctogr(xc(:,ix+1,iy+1),rlatb,rlonb)
1586  call ctogr(xc(:,ix ,iy+1),rlatc,rlonc)
1587  drlata=rlata-rlat; drlona=rlona-rlon
1588  drlatb=rlatb-rlat; drlonb=rlonb-rlon
1589  drlatc=rlatc-rlat; drlonc=rlonc-rlon
1590 
1591 ! If 'I' is the grid point (ix,iy), 'A' is (ix+1,iy); 'B' is (ix+1,iy+1)
1592 ! and 'C' is (ix,iy+1), then the area of the grid cell IABC is the sum of
1593 ! the areas of the triangles, IAB and IBC (the latter being the negative
1594 ! of the signed area of triangle, ICB):
1595  garea(ix,iy)=sarea(rlat, drlata,drlona, drlatb,drlonb) &
1596  -sarea(rlat, drlatc,drlonc, drlatb,drlonb)
1597  enddo
1598 enddo
1599 end subroutine hgrid_ak_rc
1600 
1627 subroutine hgrid_ak_dd(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, & ! [hgrid_ak_dd]
1628  delx,dely, gdlat,gdlon,garea, ff)
1629 implicit none
1630 integer(spi), intent(in ):: lx,ly,nx,ny
1631 real(dp), intent(in ):: a,k,pdlat,pdlon,&
1632  pdazi,delx,dely
1633 real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: gdlat,gdlon
1634 real(dp),dimension(lx:lx+nx-1,ly:ly+ny-1),intent(out):: garea
1635 logical, intent(out):: ff
1636 real(dp):: plat,plon,pazi
1637 plat=pdlat*dtor ! Convert these angles from degrees to radians
1638 plon=pdlon*dtor !
1639 pazi=pdazi*dtor !
1640 call hgrid_ak_rr(lx,ly,nx,ny,a,k,plat,plon,pazi, &
1641  delx,dely, gdlat,gdlon,garea, ff)
1642 if(ff)return
1643 gdlat=gdlat*rtod ! Convert these angles from radians to degrees
1644 gdlon=gdlon*rtod !
1645 end subroutine hgrid_ak_dd
1646 
1670 subroutine hgrid_ak_dd_c(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, &! [hgrid_ak_dd]
1671  delx,dely, gdlat,gdlon,garea,dx,dy,dangle_dx,dangle_dy, ff)
1672 implicit none
1673 integer(spi), intent(in ):: lx,ly,nx,ny
1674 real(dp), intent(in ):: a,k,pdlat,pdlon,&
1675  pdazi,delx,dely
1676 real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: gdlat,gdlon
1677 real(dp),dimension(lx:lx+nx-1,ly:ly+ny-1),intent(out):: garea
1678 real(dp),dimension(lx:lx+nx-1,ly:ly+ny ),intent(out):: dx
1679 real(dp),dimension(lx:lx+nx ,ly:ly+ny-1),intent(out):: dy
1680 real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: dangle_dx,dangle_dy
1681 logical, intent(out):: ff
1682 real(dp):: plat,plon,pazi
1683 plat=pdlat*dtor ! Convert these angles from degrees to radians
1684 plon=pdlon*dtor !
1685 pazi=pdazi*dtor !
1686 call hgrid_ak_rr_c(lx,ly,nx,ny,a,k,plat,plon,pazi, &
1687  delx,dely, gdlat,gdlon,garea,dx,dy,dangle_dx,dangle_dy, ff)
1688 if(ff)return
1689 gdlat =gdlat *rtod ! Convert these angles from radians to degrees
1690 gdlon =gdlon *rtod !
1691 dangle_dx=dangle_dx*rtod !
1692 dangle_dy=dangle_dy*rtod !
1693 end subroutine hgrid_ak_dd_c
1694 
1721 subroutine hgrid_ak_dc(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, & ! [hgrid_ak_dc]
1722  delx,dely, xc,xcd,garea, ff)
1723 implicit none
1724 integer(spi),intent(in ):: lx,ly,nx,ny
1725 real(dp), intent(in ):: a,k,pdlat,pdlon,pdazi,delx,dely
1726 real(dp),dimension(3, lx:lx+nx ,ly:ly+ny ),intent(out):: xc
1727 real(dp),dimension(3,2,lx:lx+nx ,ly:ly+ny ),intent(out):: xcd
1728 real(dp),dimension( lx:lx+nx-1,ly:ly+ny-1),intent(out):: garea
1729 logical, intent(out):: ff
1730 real(dp):: plat,plon,pazi
1731 plat=pdlat*dtor
1732 plon=pdlon*dtor
1733 pazi=pdazi*dtor
1734 call hgrid_ak_rc(lx,ly,nx,ny,a,k,plat,plon,pazi, &
1735  delx,dely, xc,xcd,garea, ff)
1736 end subroutine hgrid_ak_dc
1737 
1763 subroutine hgrid_ak(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak]
1764  re,delxre,delyre, glat,glon,garea, ff)
1765 implicit none
1766 integer(spi), intent(in ):: lx,ly,nx,ny
1767 real(dp), intent(in ):: a,k,plat,plon,pazi, &
1768  re,delxre,delyre
1769 real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: glat,glon
1770 real(dp),dimension(lx:lx+nx-1,ly:ly+ny-1),intent(out):: garea
1771 logical, intent(out):: ff
1772 real(dp):: delx,dely,rere
1773 delx=delxre/re ! <- set nondimensional grid step delx
1774 dely=delyre/re ! <- set nondimensional grid step dely
1775 call hgrid_ak_rr(lx,ly,nx,ny,a,k,plat,plon,pazi, &
1776  delx,dely, glat,glon,garea, ff)
1777 if(ff)return
1778 rere=re*re
1779 garea=garea*rere ! <- Convert from steradians to physical area units.
1780 end subroutine hgrid_ak
1781 
1815 subroutine hgrid_ak_c(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak]
1816  re,delxre,delyre, glat,glon,garea,dx,dy,dangle_dx,dangle_dy, ff)
1817 implicit none
1818 integer(spi), intent(in ):: lx,ly,nx,ny
1819 real(dp), intent(in ):: a,k,plat,plon,pazi, &
1820  re,delxre,delyre
1821 real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: glat,glon
1822 real(dp),dimension(lx:lx+nx-1,ly:ly+ny-1),intent(out):: garea
1823 real(dp),dimension(lx:lx+nx-1,ly:ly+ny ),intent(out):: dx
1824 real(dp),dimension(lx:lx+nx ,ly:ly+ny-1),intent(out):: dy
1825 real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: dangle_dx,dangle_dy
1826 logical, intent(out):: ff
1827 real(dp):: delx,dely,rere
1828 delx=delxre/re ! <- set nondimensional grid step delx
1829 dely=delyre/re ! <- set nondimensional grid step dely
1830 call hgrid_ak_rr_c(lx,ly,nx,ny,a,k,plat,plon,pazi, &
1831  delx,dely, glat,glon,garea,dx,dy,dangle_dx,dangle_dy, ff)
1832 if(ff)return
1833 rere=re*re
1834 garea=garea*rere ! <- Convert from steradians to physical area units.
1835 dx=dx*re ! <- Convert from nondimensional to physical length units.
1836 dy=dy*re ! <-
1837 dangle_dx=dangle_dx*rtod ! <-Convert from radians to degrees
1838 dangle_dy=dangle_dy*rtod ! <-Convert from radians to degrees
1839 end subroutine hgrid_ak_c
1840 
1850 subroutine gaulegh(m,x,w)! [gaulegh]
1851 implicit none
1852 integer(spi), intent(IN ):: m ! <- number of nodes in half-interval
1853 real(dp),dimension(m),intent(OUT):: x,w ! <- nodes and weights
1854 integer(spi),parameter:: nit=8
1855 real(dp), parameter:: eps=3.e-14_dp
1856 integer(spi) :: i,ic,j,jm,it,m2,m4p,m4p3
1857 real(dp) :: z,zzm,p1,p2,p3,pp,z1
1858 m2=m*2; m4p=m*4+1; m4p3=m4p+2
1859 do i=1,m; ic=m4p3-4*i
1860  z=cos(pih*ic/m4p); zzm=z*z-u1
1861  do it=1,nit
1862  p1=u1; p2=u0
1863  do j=1,m2; jm=j-1; p3=p2; p2=p1; p1=((j+jm)*z*p2-jm*p3)/j; enddo
1864  pp=m2*(z*p1-p2)/zzm; z1=z; z=z1-p1/pp; zzm=z*z-u1
1865  if(abs(z-z1) <= eps)exit
1866  enddo
1867  x(i)=z; w(i)=-u2/(zzm*pp*pp)
1868 enddo
1869 end subroutine gaulegh
1870 
1886 subroutine gtoxm_ak_rr_m(A,K,plat,plon,pazi,lat,lon,xm,ff)! [gtoxm_ak_rr]
1887 use pmat5, only: grtoc
1888 implicit none
1889 real(dp), intent(in ):: a,k,plat,plon,pazi,lat,lon
1890 real(dp),dimension(2),intent(out):: xm
1891 logical, intent(out):: ff
1892 real(dp),dimension(3,3):: prot,azirot
1893 real(dp) :: clat,slat,clon,slon,cazi,sazi
1894 real(dp),dimension(3) :: xc
1895 clat=cos(plat); slat=sin(plat)
1896 clon=cos(plon); slon=sin(plon)
1897 cazi=cos(pazi); sazi=sin(pazi)
1898 
1899 azirot(:,1)=(/ cazi, sazi, u0/)
1900 azirot(:,2)=(/-sazi, cazi, u0/)
1901 azirot(:,3)=(/ u0, u0, u1/)
1902 
1903 prot(:,1)=(/ -slon, clon, u0/)
1904 prot(:,2)=(/-slat*clon, -slat*slon, clat/)
1905 prot(:,3)=(/ clat*clon, clat*slon, slat/)
1906 prot=matmul(prot,azirot)
1907 
1908 call grtoc(lat,lon,xc)
1909 xc=matmul(transpose(prot),xc)
1910 call xctoxm_ak(a,k,xc,xm,ff)
1911 end subroutine gtoxm_ak_rr_m
1912 
1930 subroutine gtoxm_ak_rr_g(A,K,plat,plon,pazi,delx,dely,lat,lon,&! [gtoxm_ak_rr]
1931  xm,ff)
1932 implicit none
1933 real(dp), intent(in ):: a,k,plat,plon,pazi,delx,dely,lat,lon
1934 real(dp),dimension(2),intent(out):: xm
1935 logical, intent(out):: ff
1936 call gtoxm_ak_rr_m(a,k,plat,plon,pazi,lat,lon,xm,ff); if(ff)return
1937 xm(1)=xm(1)/delx; xm(2)=xm(2)/dely
1938 end subroutine gtoxm_ak_rr_g
1939 
1952 subroutine gtoxm_ak_dd_m(A,K,pdlat,pdlon,pdazi,dlat,dlon,&! [gtoxm_ak_dd]
1953  xm,ff)
1954 implicit none
1955 real(dp), intent(in ):: a,k,pdlat,pdlon,pdazi,dlat,dlon
1956 real(dp),dimension(2),intent(out):: xm
1957 logical, intent(out):: ff
1958 real(dp):: plat,plon,pazi,lat,lon
1959 plat=pdlat*dtor ! Convert these angles from degrees to radians
1960 plon=pdlon*dtor !
1961 pazi=pdazi*dtor !
1962 lat=dlat*dtor
1963 lon=dlon*dtor
1964 call gtoxm_ak_rr_m(a,k,plat,plon,pazi,lat,lon,xm,ff)
1965 end subroutine gtoxm_ak_dd_m
1966 
1981 subroutine gtoxm_ak_dd_g(A,K,pdlat,pdlon,pdazi,delx,dely,&! [gtoxm_ak_dd]
1982 dlat,dlon, xm,ff)
1983 implicit none
1984 real(dp), intent(in ):: a,k,pdlat,pdlon,pdazi,delx,dely,dlat,dlon
1985 real(dp),dimension(2),intent(out):: xm
1986 logical, intent(out):: ff
1987 real(dp):: plat,plon,pazi,lat,lon
1988 plat=pdlat*dtor ! Convert these angles from degrees to radians
1989 plon=pdlon*dtor !
1990 pazi=pdazi*dtor !
1991 lat=dlat*dtor
1992 lon=dlon*dtor
1993 call gtoxm_ak_rr_g(a,k,plat,plon,pazi,delx,dely,lat,lon,xm,ff)
1994 end subroutine gtoxm_ak_dd_g
1995 
2013 subroutine xmtog_ak_rr_m(A,K,plat,plon,pazi,xm,lat,lon,ff)! [xmtog_ak_rr]
2014 use pmat5, only: ctogr
2015 implicit none
2016 real(dp), intent(in ):: a,k,plat,plon,pazi
2017 real(dp),dimension(2),intent(in ):: xm
2018 real(dp), intent(out):: lat,lon
2019 logical, intent(out):: ff
2020 real(dp),dimension(3,2):: xcd
2021 real(dp),dimension(3,3):: prot,azirot
2022 real(dp) :: clat,slat,clon,slon,cazi,sazi
2023 real(dp),dimension(3) :: xc
2024 clat=cos(plat); slat=sin(plat)
2025 clon=cos(plon); slon=sin(plon)
2026 cazi=cos(pazi); sazi=sin(pazi)
2027 
2028 azirot(:,1)=(/ cazi, sazi, u0/)
2029 azirot(:,2)=(/-sazi, cazi, u0/)
2030 azirot(:,3)=(/ u0, u0, u1/)
2031 
2032 prot(:,1)=(/ -slon, clon, u0/)
2033 prot(:,2)=(/-slat*clon, -slat*slon, clat/)
2034 prot(:,3)=(/ clat*clon, clat*slon, slat/)
2035 prot=matmul(prot,azirot)
2036 call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
2037 xc=matmul(prot,xc)
2038 call ctogr(xc,lat,lon)
2039 end subroutine xmtog_ak_rr_m
2040 
2060 subroutine xmtog_ak_rr_g(A,K,plat,plon,pazi,delx,dely,xm,&! [xmtog_ak_rr]
2061  lat,lon,ff)
2062 implicit none
2063 real(dp), intent(in ):: a,k,plat,plon,pazi,delx,dely
2064 real(dp),dimension(2),intent(in ):: xm
2065 real(dp), intent(out):: lat,lon
2066 logical, intent(out):: ff
2067 real(dp),dimension(2):: xmt
2068 xmt(1)=xm(1)*delx ! Convert from grid units to intrinsic map-space units
2069 xmt(2)=xm(2)*dely !
2070 call xmtog_ak_rr_m(a,k,plat,plon,pazi,xmt,lat,lon,ff)
2071 end subroutine xmtog_ak_rr_g
2072 
2085 subroutine xmtog_ak_dd_m(A,K,pdlat,pdlon,pdazi,xm,dlat,dlon,ff)! [xmtog_ak_dd]
2086 use pmat5, only: ctogr
2087 implicit none
2088 real(dp), intent(in ):: a,k,pdlat,pdlon,pdazi
2089 real(dp),dimension(2),intent(in ):: xm
2090 real(dp), intent(out):: dlat,dlon
2091 logical, intent(out):: ff
2092 real(dp):: plat,plon,pazi,lat,lon
2093 plat=pdlat*dtor ! Convert these angles from degrees to radians
2094 plon=pdlon*dtor !
2095 pazi=pdazi*dtor !
2096 call xmtog_ak_rr_m(a,k,plat,plon,pazi,xm,lat,lon,ff)
2097 dlat=lat*rtod
2098 dlon=lon*rtod
2099 end subroutine xmtog_ak_dd_m
2100 
2115 subroutine xmtog_ak_dd_g(A,K,pdlat,pdlon,pdazi,delx,dely,xm,&! [xmtog_ak_dd]
2116  dlat,dlon,ff)
2117 implicit none
2118 real(dp), intent(in ):: a,k,pdlat,pdlon,pdazi,delx,dely
2119 real(dp),dimension(2),intent(in ):: xm
2120 real(dp), intent(out):: dlat,dlon
2121 logical, intent(out):: ff
2122 real(dp),dimension(2):: xmt
2123 real(dp) :: plat,plon,pazi,lat,lon
2124 xmt(1)=xm(1)*delx ! Convert from grid units to intrinsic map-space units
2125 xmt(2)=xm(2)*dely !
2126 plat=pdlat*dtor ! Convert these angles from degrees to radians
2127 plon=pdlon*dtor !
2128 pazi=pdazi*dtor !
2129 call xmtog_ak_rr_m(a,k,plat,plon,pazi,xmt,lat,lon,ff)
2130 dlat=lat*rtod
2131 dlon=lon*rtod
2132 end subroutine xmtog_ak_dd_g
2133 
2134 end module pesg
2135 
pesg::xstoxc::xstoxc1
subroutine xstoxc1(xs, xc, xcd, xcdd)
Standard transformation from polar stereographic map coordinates, xs, to cartesian, xc, assuming the projection axis is polar.
Definition: pesg.f90:106
pesg::xmtoxc_ak::xmtoxc_vak1
subroutine xmtoxc_vak1(ak, xm, xc, xcd, xc1, xcd1, ff)
Like xmtoxc_vak, _ak, but also return derivatives wrt ak.
Definition: pesg.f90:399
pesg::xmtoxc_ak::xmtoxc_vak
subroutine xmtoxc_vak(ak, xm, xc, xcd, ff)
Assuming the vector AK parameterization of the Extended Schmidt-transformed Gnomonic (ESG) mapping wi...
Definition: pesg.f90:380
pesg::xmtog_ak_rr::xmtog_ak_rr_m
subroutine xmtog_ak_rr_m(A, K, plat, plon, pazi, xm, lat, lon, ff)
Given the ESG map specified by parameters (A,K) and geographical center and orientation, plat,plon,pazi (radians), and a position, in map-space coordinates given by the 2-vector, xm, return the geographical coordinates, lat and lon (radians).
Definition: pesg.f90:2013
pmat5
Utility routines for orienting the globe and basic geographical mappings.
Definition: pmat5.f90:27
pesg::xctoxm_ak
Definition: pesg.f90:33
pesg::xmtog_ak_dd
Definition: pesg.f90:55
pmat4::outer_product
Definition: pmat4.f90:67
pesg::hgrid_ak
Definition: pesg.f90:48
pesg::get_meanq::get_meanqd
subroutine get_meanqd(ngh, lam, xg, wg, ak, ma, q, qdak, qdma, ga, gadak, gadma, ff)
For a parameter vector, ak and a map-space domain-parameter vector, ma, return the lambda-parameteriz...
Definition: pesg.f90:601
pesg::bestesg_geo
Definition: pesg.f90:42
pesg::xmtog_ak_rr::xmtog_ak_rr_g
subroutine xmtog_ak_rr_g(A, K, plat, plon, pazi, delx, dely, xm, lat, lon, ff)
For an ESG map with parameters, (A,K), and geographical orientation, given by plon,plat,pazi (radians), and given a point in grid-space units as the 2-vector, xm, return the geographical coordinates, lat, lon, (radians) of this point.
Definition: pesg.f90:2060
pesg::hgrid_ak_rr
Definition: pesg.f90:44
pesg::zmtozt::zmtozt1
subroutine zmtozt1(a, zm, zt, ztd, zt1, ztd1, ff)
Like zmtozt, but also, get the derivative with respect to a, zt1 of zt, and ztd1 of ztd...
Definition: pesg.f90:347
pesg::get_qofv::get_qsofvs
subroutine get_qsofvs(n, lam, v1s, v2s, v3s, v4s, qs)
General util to convert value.
Definition: pesg.f90:769
pesg::get_qx
Definition: pesg.f90:37
pesg::gtoxm_ak_dd::gtoxm_ak_dd_g
subroutine gtoxm_ak_dd_g(A, K, pdlat, pdlon, pdazi, delx, dely, dlat, dlon, xm, ff)
Like gtoxm_ak_rr_g, except lat, lon, azimuth, are expressed in degrees.
Definition: pesg.f90:1981
pmat4::cross_product
Definition: pmat4.f90:64
pesg::hgrid_ak_dd::hgrid_ak_dd_c
subroutine hgrid_ak_dd_c(lx, ly, nx, ny, a, k, pdlat, pdlon, pdazi, delx, dely, gdlat, gdlon, garea, dx, dy, dangle_dx, dangle_dy, ff)
Like hgrid_ak_rr_c, except all the angle arguments (but not delx,dely) are in degrees instead of radi...
Definition: pesg.f90:1670
pesg::gtoxm_ak_rr::gtoxm_ak_rr_m
subroutine gtoxm_ak_rr_m(A, K, plat, plon, pazi, lat, lon, xm, ff)
Given the map specification (angles in radians), the grid spacing (in map-space units) and the sample...
Definition: pesg.f90:1886
pmat4
Module for handy vector and matrix operations in Euclidean geometry.
Definition: pmat4.f90:45
pesg::gtoxm_ak_dd::gtoxm_ak_dd_m
subroutine gtoxm_ak_dd_m(A, K, pdlat, pdlon, pdazi, dlat, dlon, xm, ff)
Like gtoxm_ak_rr_m, except lat, lon, azimuth, are expressed in degrees.
Definition: pesg.f90:1952
pfun
This module is for evaluating several useful real-valued functions that are not always available in F...
Definition: pfun.f90:11
pesg::get_meanq
Definition: pesg.f90:39
pesg::xmtog_ak_rr
Definition: pesg.f90:53
pesg
Suite of routines to perform the 2-parameter family of Extended Schmidt Gnomonic (ESG) regional grid ...
Definition: pesg.f90:15
pkind
Standard integer, real, and complex single and double precision kinds.
Definition: pkind.f90:7
pesg::hgrid_ak_rc
Definition: pesg.f90:45
pesg::gtoxm_ak_dd
Definition: pesg.f90:51
pesg::bestesg_map
Definition: pesg.f90:43
pesg::xmtog_ak_dd::xmtog_ak_dd_m
subroutine xmtog_ak_dd_m(A, K, pdlat, pdlon, pdazi, xm, dlat, dlon, ff)
Like xmtog_ak_rr_m, except lat, lon, azimuth, are expressed in degrees.
Definition: pesg.f90:2085
pesg::hgrid_ak::hgrid_ak_c
subroutine hgrid_ak_c(lx, ly, nx, ny, a, k, plat, plon, pazi, re, delxre, delyre, glat, glon, garea, dx, dy, dangle_dx, dangle_dy, ff)
Like hgrid_ak_rr_c except the argument list includes the earth radius, re, and this is used to expres...
Definition: pesg.f90:1815
pmat4::triple_product
Definition: pmat4.f90:70
psym2
A suite of routines to perform the eigen-decomposition of symmetric 2*2 matrices and to deliver basic...
Definition: psym2.f90:11
pietc
Some of the commonly used constants (pi etc) mainly for double-precision subroutines.
Definition: pietc.f90:14
pesg::hgrid_ak_dc
Definition: pesg.f90:47
pesg::get_edges
Definition: pesg.f90:36
pmat
Definition: pmat.f90:20
pesg::get_meanq::get_meanqs
subroutine get_meanqs(n, ngh, lam, xg, wg, aks, mas, qs, ff)
Like getmeanqd, except for n different values, aks, of ak and n different values, mas of ma...
Definition: pesg.f90:680
pesg::guessak_map
Definition: pesg.f90:40
pesg::xmtoxc_ak
Definition: pesg.f90:34
pesg::gtoxm_ak_rr::gtoxm_ak_rr_g
subroutine gtoxm_ak_rr_g(A, K, plat, plon, pazi, delx, dely, lat, lon, xm, ff)
Given the map specification (angles in radians), the grid spacing (in map-space units) and the sample...
Definition: pesg.f90:1930
pesg::get_qofv
Definition: pesg.f90:38
pesg::xmtoxt::xmtoxt1
subroutine xmtoxt1(a, xm, xt, xtd, xt1, xtd1, ff)
Like zmtozt1, but for 2-vector xm and xt, and 2*2 diagonal Jacobian xtd Also, the derivatives...
Definition: pesg.f90:275
pesg::hgrid_ak_rr::hgrid_ak_rr_c
subroutine hgrid_ak_rr_c(lx, ly, nx, ny, a, k, plat, plon, pazi, delx, dely, glat, glon, garea, dx, dy, angle_dx, angle_dy, ff)
Use a and k as the parameters of an extended Schmidt-transformed gnomonic (ESG) mapping centered at (...
Definition: pesg.f90:1335
pesg::guessak_geo
Definition: pesg.f90:41
pesg::hgrid_ak_dd
Definition: pesg.f90:46
pesg::xttoxs::xttoxs1
subroutine xttoxs1(k, xt, xs, xsd, xsdd, xs1, xsd1, ff)
Like xttoxs, but also, return the derivatives, wrt K, of xs and xsd.
Definition: pesg.f90:192
pesg::xmtog_ak_dd::xmtog_ak_dd_g
subroutine xmtog_ak_dd_g(A, K, pdlat, pdlon, pdazi, delx, dely, xm, dlat, dlon, ff)
Like xmtog_ak_rr_g, except lat, lon, azimuth, are expressed in degrees.
Definition: pesg.f90:2115
pesg::get_qofv::get_qofvd
subroutine get_qofvd(lam, v2, v3, v1d, v2d, v3d, v4d, qd)
Like get_qofv, but for (only) the 2-vector derivatives of Q.
Definition: pesg.f90:749
pesg::get_qx::get_qxd
subroutine get_qxd(j0, j0d, v1, v2, v3, v4, v1d, v2d, v3d, v4d)
From a jacobian matrix, j0, and its derivative, j0d, get a sufficient set of v.
Definition: pesg.f90:540
pesg::gtoxm_ak_rr
Definition: pesg.f90:49