grid_tools  1.13.0
 All Data Structures Files Functions Variables Pages
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,mat22
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 mat22=xsdd(:,:,1)*xt1(1)+xsdd(:,:,2)*xt1(2)
416 xsd1(:,:,1)=matmul(xsd,xtd1)+matmul(mat22,xtd)
417 xsd1(:,:,2)=matmul(xsdk,xtd)
418 xsd=matmul(xsd,xtd)
419 call xstoxc(xs,xc,xcd,xcdd)
420 xc1=matmul(xcd,xs1)
421 do i=1,3; xcd1(i,:,:)=matmul(transpose(xsd),matmul(xcdd(i,:,:),xs1)); enddo
422 do i=1,2; xcd1(:,:,i)=xcd1(:,:,i)+matmul(xcd,xsd1(:,:,i)); enddo
423 xcd=matmul(xcd,xsd)
424 end subroutine xmtoxc_vak1
425 
439 subroutine xmtoxc_ak(a,k,xm,xc,xcd,ff)! [xmtoxc_ak]
440 implicit none
441 real(dp), intent(in ):: a,k
442 real(dp),dimension(2), intent(in ):: xm
443 real(dp),dimension(3), intent(out):: xc
444 real(dp),dimension(3,2),intent(out):: xcd
445 logical, intent(out):: ff
446 real(dp),dimension(2,2):: xtd,xsd
447 real(dp),dimension(2) :: xt,xs
448 call xmtoxt(a,xm,xt,xtd,ff); if(ff)return
449 call xttoxs(k,xt,xs,xsd,ff); if(ff)return
450 xsd=matmul(xsd,xtd)
451 call xstoxc(xs,xc,xcd)
452 xcd=matmul(xcd,xsd)
453 end subroutine xmtoxc_ak
454 
465 subroutine xctoxm_ak(a,k,xc,xm,ff)! [xctoxm_ak]
466 implicit none
467 real(dp), intent(in ):: a,k
468 real(dp),dimension(3),intent(in ):: xc
469 real(dp),dimension(2),intent(out):: xm
470 logical, intent(out):: ff
471 real(dp),dimension(2):: xs,xt
472 ff=f
473 call xctoxs(xc,xs)
474 call xstoxt(k,xs,xt,ff); if(ff)return
475 call xttoxm(a,xt,xm,ff)
476 end subroutine xctoxm_ak
477 
488 subroutine get_edges(arcx,arcy,edgex,edgey)! [get_edges]
489 implicit none
490 real(dp), intent(in ):: arcx,arcy
491 real(dp),dimension(3),intent(out):: edgex,edgey
492 real(dp):: cx,sx,cy,sy,rarcx,rarcy
493 rarcx=arcx*dtor; rarcy=arcy*dtor
494 cx=cos(rarcx); sx=sin(rarcx)
495 cy=cos(rarcy); sy=sin(rarcy)
496 edgex=(/sx,u0,cx/); edgey=(/u0,sy,cy/)
497 end subroutine get_edges
498 
511 subroutine get_qx(j0, v1,v2,v3,v4)! [get_qx]
512 use psym2, only: logsym2
513 implicit none
514 real(dp),dimension(3,2),intent(in ):: j0
515 real(dp), intent(out):: v1,v2,v3,v4
516 real(dp),dimension(2,2):: el,g
517 g=matmul(transpose(j0),j0)
518 call logsym2(g,el); el=el*o2
519 v1=el(1,1)**2+u2*el(1,2)**2+el(2,2)**2
520 v2=el(1,1)
521 v3=el(2,2)
522 v4=(el(1,1)+el(2,2))**2
523 end subroutine get_qx
524 
541 subroutine get_qxd(j0,j0d, v1,v2,v3,v4,v1d,v2d,v3d,v4d)! [get_qx]
542 use psym2, only: logsym2
543 implicit none
544 real(dp),dimension(3,2), intent(in ):: j0
545 real(dp),dimension(3,2,2),intent(in ):: j0d
546 real(dp), intent(out):: v1,v2,v3,v4
547 real(dp),dimension(2), intent(out):: v1d,v2d,v3d,v4d
548 real(dp),dimension(2,2,2,2):: deldg
549 real(dp),dimension(2,2,2) :: eld,gd
550 real(dp),dimension(2,2) :: el,g
551 integer(spi) :: i,j,k
552 g=matmul(transpose(j0),j0)
553 do i=1,2
554  gd(:,:,i)=matmul(transpose(j0d(:,:,i)),j0)+matmul(transpose(j0),j0d(:,:,i))
555 enddo
556 call logsym2(g,el,deldg); el=el*o2; deldg=deldg*o2
557 eld=u0
558 do i=1,2; do j=1,2; do k=1,2
559  eld(:,:,k)=eld(:,:,k)+deldg(:,:,i,j)*gd(i,j,k)
560 enddo ; enddo ; enddo
561 v1=el(1,1)**2+u2*el(1,2)**2+el(2,2)**2
562 v2=el(1,1)
563 v3=el(2,2)
564 v4=(el(1,1)+el(2,2))**2
565 v1d=u2*(el(1,1)*eld(1,1,:)+u2*el(1,2)*eld(1,2,:)+el(2,2)*eld(2,2,:))
566 v2d=eld(1,1,:)
567 v3d=eld(2,2,:)
568 v4d=u2*(el(1,1)+el(2,2))*(eld(1,1,:)+eld(2,2,:))
569 end subroutine get_qxd
570 
602 subroutine get_meanqd(ngh,lam,xg,wg,ak,ma, q,qdak,qdma, & ! [get_meanq]
603  ga,gadak,gadma, ff)
604 implicit none
605 integer(spi), intent(in ):: ngh
606 real(dp), intent(in ):: lam
607 real(dp),dimension(ngh),intent(in ):: xg,wg
608 real(dp),dimension(2) ,intent(in ):: ak,ma
609 real(dp), intent(out):: q
610 real(dp),dimension(2), intent(out):: qdak,qdma
611 real(dp),dimension(2), intent(out):: ga
612 real(dp),dimension(2,2),intent(out):: gadak,gadma
613 logical, intent(out):: ff
614 real(dp),dimension(3,2,2):: xcd1
615 real(dp),dimension(3,2) :: xcd,xc1
616 real(dp),dimension(3) :: xc
617 real(dp),dimension(2) :: xm, v1dxy,v2dxy,v3dxy,v4dxy, &
618  v1dl,v2dl,v3dl,v4dl, v1d,v2d,v3d,v4d
619 real(dp) :: wx,wy, &
620  v1xy,v2xy,v3xy,v4xy, v1l,v2l,v3l,v4l, v1,v2,v3,v4
621 integer(spi) :: i,ic,ix,iy
622 v1 =u0; v2 =u0; v3 =u0; v4 =u0
623 v1d=u0; v2d=u0; v3d=u0; v4d=u0
624 do iy=1,ngh
625  wy=wg(iy)
626  xm(2)=ma(2)*xg(iy)
627  v1l =u0; v2l =u0; v3l =u0; v4l =u0
628  v1dl=u0; v2dl=u0; v3dl=u0; v4dl=u0
629  do ix=1,ngh
630  wx=wg(ix)
631  xm(1)=ma(1)*xg(ix)
632  call xmtoxc_ak(ak,xm,xc,xcd,xc1,xcd1,ff); if(ff)return
633  call get_qx(xcd,xcd1, v1xy,v2xy,v3xy,v4xy, v1dxy,v2dxy,v3dxy,v4dxy)
634  v1l =v1l +wx*v1xy; v2l =v2l +wx*v2xy
635  v3l =v3l +wx*v3xy; v4l =v4l +wx*v4xy
636  v1dl=v1dl+wx*v1dxy;v2dl=v2dl+wx*v2dxy
637  v3dl=v3dl+wx*v3dxy;v4dl=v4dl+wx*v4dxy
638  enddo
639  v1 =v1 +wy*v1l; v2 =v2 +wy*v2l; v3 =v3 +wy*v3l; v4 =v4 +wy*v4l
640  v1d=v1d+wy*v1dl; v2d=v2d+wy*v2dl; v3d=v3d+wy*v3dl; v4d=v4d+wy*v4dl
641 enddo
642 call get_qofv(lam,v1,v2,v3,v4, q)! <- Q(lam) based on the v1,v2,v3,v4
643 call get_qofv(lam,v2,v3, v1d,v2d,v3d,v4d, qdak)! <- Derivative of Q wrt ak
644 ! Derivatives of ga wrt ak, and of q and ga wrt ma:
645 gadma=u0! <- needed because only diagonal elements are filled
646 do i=1,2
647  ic=3-i
648  xm=0; xm(i)=ma(i)
649  call xmtoxc_ak(ak,xm,xc,xcd,xc1,xcd1,ff); if(ff)return
650  ga(i)=atan2(xc(i),xc(3))*rtod
651  gadak(i,:)=(xc(3)*xc1(i,:)-xc(i)*xc1(3,:))*rtod
652  gadma(i,i)=(xc(3)*xcd(i,i)-xc(i)*xcd(3,i))*rtod
653 
654  v1l=u0; v2l=u0; v3l=u0; v4l=u0
655  do iy=1,ngh
656  wy=wg(iy)
657  xm(ic)=ma(ic)*xg(iy)
658  call xmtoxc_ak(ak,xm,xc,xcd,ff); if(ff)return
659  call get_qx(xcd, v1xy,v2xy,v3xy,v4xy)
660  v1l=v1l+wy*v1xy; v2l=v2l+wy*v2xy; v3l=v3l+wy*v3xy; v4l=v4l+wy*v4xy
661  enddo
662  v1d(i)=(v1l-v1)/ma(i); v2d(i)=(v2l-v2)/ma(i)
663  v3d(i)=(v3l-v3)/ma(i); v4d(i)=(v4l-v4)/ma(i)
664 enddo
665 call get_qofv(lam,v2,v3, v1d,v2d,v3d,v4d, qdma)
666 end subroutine get_meanqd
667 
681 subroutine get_meanqs(n,ngh,lam,xg,wg,aks,mas, qs,ff)! [get_meanq]
682 implicit none
683 integer(spi), intent(in ):: n,ngh
684 real(dp),dimension(ngh),intent(in ):: xg,wg
685 real(dp), intent(in ):: lam
686 real(dp),dimension(2,n),intent(in ):: aks,mas
687 real(dp),dimension(n), intent(out):: qs
688 logical, intent(out):: ff
689 real(dp),dimension(n) :: v1s,v2s,v3s,v4s
690 real(dp),dimension(n) :: v1sl,v2sl,v3sl,v4sl
691 real(dp),dimension(3,2):: xcd
692 real(dp),dimension(3) :: xc
693 real(dp),dimension(2) :: xm,xgs
694 real(dp) :: wx,wy, v1xy,v2xy,v3xy,v4xy
695 integer(spi) :: i,ix,iy
696 v1s=u0; v2s=u0; v3s=u0; v4s=u0
697 do iy=1,ngh
698  wy=wg(iy)
699  v1sl=u0; v2sl=u0; v3sl=u0; v4sl=u0
700  do ix=1,ngh
701  wx=wg(ix)
702  xgs=(/xg(ix),xg(iy)/)
703  do i=1,n
704  xm=mas(:,i)*xgs
705  call xmtoxc_ak(aks(:,i),xm,xc,xcd,ff); if(ff)return
706  call get_qx(xcd,v1xy,v2xy,v3xy,v4xy)
707  v1sl(i)=v1sl(i)+wx*v1xy; v2sl(i)=v2sl(i)+wx*v2xy
708  v3sl(i)=v3sl(i)+wx*v3xy; v4sl(i)=v4sl(i)+wx*v4xy
709  enddo
710  enddo
711  v1s=v1s+wy*v1sl; v2s=v2s+wy*v2sl; v3s=v3s+wy*v3sl; v4s=v4s+wy*v4sl
712 enddo
713 call get_qofv(n,lam,v1s,v2s,v3s,v4s, qs)
714 end subroutine get_meanqs
715 
728 subroutine get_qofv(lam,v1,v2,v3,v4, q)! [get_qofv]
729 implicit none
730 real(dp),intent(in ):: lam,v1,v2,v3,v4
731 real(dp),intent(out):: q
732 real(dp):: lamc
733 lamc=u1-lam
734 q=lamc*(v1-(v2**2+v3**2)) +lam*(v4 -(v2+v3)**2)
735 end subroutine get_qofv
736 
750 subroutine get_qofvd(lam, v2,v3, v1d,v2d,v3d,v4d, qd)! [get_qofv]
751 implicit none
752 real(dp), intent(in ):: lam,v2,v3
753 real(dp),dimension(2),intent(in ):: v1d,v2d,v3d,v4d
754 real(dp),dimension(2),intent(out):: qd
755 real(dp):: lamc
756 lamc=u1-lam
757 qd=lamc*(v1d-u2*(v2d*v2+v3d*v3))+lam*(v4d-u2*(v2d+v3d)*(v2+v3))
758 end subroutine get_qofvd
759 
770 subroutine get_qsofvs(n,lam,v1s,v2s,v3s,v4s, qs)! [get_qofv]
771 implicit none
772 integer(spi), intent(in ):: n
773 real(dp), intent(in ):: lam
774 real(dp),dimension(n),intent(in ):: v1s,v2s,v3s,v4s
775 real(dp),dimension(n),intent(out):: qs
776 real(dp):: lamc
777 lamc=u1-lam
778 qs=lamc*(v1s-(v2s**2+v3s**2)) +lam*(v4s -(v2s+v3s)**2)
779 end subroutine get_qsofvs
780 
790 subroutine guessak_map(asp,tmarcx,ak)! [guessak_map]
791 implicit none
792 real(dp), intent(in ):: asp,tmarcx
793 real(dp),dimension(2),intent(out):: ak
794 real(dp):: gmarcx
795 gmarcx=tmarcx*rtod
796 call guessak_geo(asp,gmarcx,ak)
797 end subroutine guessak_map
798 
808 subroutine guessak_geo(asp,arc,ak)! [guessak_geo]
809 implicit none
810 real(dp), intent(in ):: asp,arc
811 real(dp),dimension(2),intent(out):: ak
812 integer(spi),parameter :: narc=11,nasp=10! <- Table index bounds
813 real(dp), parameter :: eps=1.e-7_dp,darc=10._dp+eps,dasp=.1_dp+eps
814 real(dp),dimension(nasp,0:narc):: adarc,kdarc
815 real(dp) :: sasp,sarc,wx0,wx1,wa0,wa1
816 integer(spi) :: iasp0,iasp1,iarc0,iarc1
817 !------------------------
818 ! Tables of approximate A (adarc) and K (kdarc), valid for lam=0.8, for aspect ratio,
819 ! asp, at .1, .2, .3, .4, .5, .6, .7, .8, 1. and major semi-arc, arcx, at values,
820 ! 0 (nominally), 10., ..., 100. degrees (where the nominal "0" angle is actually
821 ! from a computation at 3 degrees, since zero would not make sense).
822 ! The 100 degree rows are repeated as the 110 degree entries deliberately to pad the
823 ! table into the partly forbidden parameter space to allow skinny domains of up to
824 ! 110 degree semi-length to be validly endowed with optimal A and K:
825 data adarc/ &
826 -.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,&
827 -.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,&
828 -.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,&
829 -.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,&
830 -.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,&
831 -.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,&
832 -.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,&
833 -.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,&
834 -.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,&
835 -.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,&
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 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/
838 
839 data kdarc/ &
840 -.947_dp,-.818_dp,-.668_dp,-.535_dp,-.433_dp,-.361_dp,-.313_dp,-.284_dp,-.269_dp,-.264_dp,&
841 -.946_dp,-.816_dp,-.665_dp,-.533_dp,-.431_dp,-.359_dp,-.311_dp,-.282_dp,-.267_dp,-.262_dp,&
842 -.942_dp,-.806_dp,-.655_dp,-.524_dp,-.424_dp,-.353_dp,-.305_dp,-.276_dp,-.261_dp,-.255_dp,&
843 -.932_dp,-.789_dp,-.637_dp,-.509_dp,-.412_dp,-.343_dp,-.296_dp,-.266_dp,-.250_dp,-.244_dp,&
844 -.909_dp,-.759_dp,-.609_dp,-.486_dp,-.394_dp,-.328_dp,-.283_dp,-.254_dp,-.237_dp,-.230_dp,&
845 -.863_dp,-.711_dp,-.569_dp,-.456_dp,-.372_dp,-.310_dp,-.266_dp,-.238_dp,-.220_dp,-.212_dp,&
846 -.779_dp,-.642_dp,-.518_dp,-.419_dp,-.343_dp,-.287_dp,-.247_dp,-.220_dp,-.202_dp,-.192_dp,&
847 -.661_dp,-.556_dp,-.456_dp,-.374_dp,-.310_dp,-.262_dp,-.226_dp,-.200_dp,-.182_dp,-.171_dp,&
848 -.533_dp,-.462_dp,-.388_dp,-.325_dp,-.274_dp,-.234_dp,-.203_dp,-.179_dp,-.161_dp,-.150_dp,&
849 -.418_dp,-.373_dp,-.322_dp,-.275_dp,-.236_dp,-.204_dp,-.178_dp,-.156_dp,-.139_dp,-.127_dp,&
850 -.324_dp,-.296_dp,-.261_dp,-.229_dp,-.200_dp,-.174_dp,-.152_dp,-.133_dp,-.117_dp,-.104_dp,&
851 -.324_dp,-.296_dp,-.261_dp,-.229_dp,-.200_dp,-.174_dp,-.152_dp,-.133_dp,-.117_dp,-.104_dp/
852 !=============================================================================
853 sasp=asp/dasp
854 iasp0=floor(sasp); wa1=sasp-iasp0
855 iasp1=iasp0+1; wa0=iasp1-sasp
856 sarc=arc/darc
857 iarc0=floor(sarc); wx1=sarc-iarc0
858 iarc1=iarc0+1; wx0=iarc1-sarc
859 if(iasp0<1 .or. iasp1>nasp)stop 'Guessak_geo; Aspect ratio out of range'
860 if(iarc0<0 .or. iarc1>narc)stop 'Guessak_geo; Major semi-arc is out of range'
861 
862 ! Bilinearly interpolate A and K from crude table into a 2-vector:
863 ak=(/wx0*(wa0*adarc(iasp0,iarc0)+wa1*adarc(iasp1,iarc0))+ &
864  wx1*(wa0*adarc(iasp0,iarc1)+wa1*adarc(iasp1,iarc1)), &
865  wx0*(wa0*kdarc(iasp0,iarc0)+wa1*kdarc(iasp1,iarc0))+ &
866  wx1*(wa0*kdarc(iasp0,iarc1)+wa1*kdarc(iasp1,iarc1))/)
867 end subroutine guessak_geo
868 
902 subroutine bestesg_geo(lam,garcx,garcy, a,k,marcx,marcy,q,ff)! [bestesg_geo]
903 use pietc, only: u5,o5,s18,s36,s54,s72,ms18,ms36,ms54,ms72
904 use pmat, only: inv
905 use psym2, only: chol2
906 implicit none
907 real(dp),intent(in ):: lam,garcx,garcy
908 real(dp),intent(out):: a,k,marcx,marcy,q
909 logical ,intent(out):: ff
910 integer(spi),parameter :: nit=200,mit=20,ngh=25
911 real(dp) ,parameter :: u2o5=u2*o5,&
912  f18=u2o5*s18,f36=u2o5*s36,f54=u2o5*s54,&
913  f72=u2o5*s72,mf18=-f18,mf36=-f36,mf54=-f54,&
914  mf72=-f72,& !<- (Fourier transform coefficients)
915  r=0.001_dp,rr=r*r,dang=pi2*o5,crit=1.e-14_dp
916 real(dp),dimension(ngh) :: wg,xg
917 real(dp),dimension(0:4,0:4):: em5 ! <- Fourier matrix for 5 points
918 real(dp),dimension(0:4) :: qs
919 real(dp),dimension(2,0:4) :: aks,mas
920 real(dp),dimension(2,2) :: basis0,basis,hess,el,gadak,gadma,madga,madak
921 real(dp),dimension(2) :: ak,dak,dma,vec2,grad,qdak,qdma,ga,ma,gat
922 real(dp) :: s,tgarcx,tgarcy,asp,ang
923 integer(spi) :: i,it
924 logical :: flip
925 data em5/o5,u2o5, u0,u2o5, u0,& ! <-The Fourier matrix for 5 points. Applied
926  o5, f18, f72,mf54, f36,& ! to the five 72-degree spaced values in a
927  o5,mf54, f36, f18,mf72,& ! column-vector, the product vector has the
928  o5,mf54,mf36, f18, f72,& ! components, wave-0, cos and sin wave-1,
929  o5, f18,mf72,mf54,mf36/ ! cos and sin wave-2.
930 ! First guess upper-triangular basis regularizing the samples used to
931 ! estimate the Hessian of q by finite differencing:
932 data basis0/0.1_dp,u0, 0.3_dp,0.3_dp/
933 ff=lam<u0 .or. lam>=u1
934 if(ff)then; print'("In bestesg_geo; lam out of range")';return; endif
935 ff= garcx<=u0 .or. garcy<=u0
936 if(ff)then
937  print'("In bestesg_geo; a nonpositive domain parameter, garcx or garcy")'
938  return
939 endif
940 call gaulegh(ngh,xg,wg)! <- Prepare Gauss-Legendre nodes xg and weights wg
941 flip=garcy>garcx
942 if(flip)then; tgarcx=garcy; tgarcy=garcx! <- Switch
943 else ; tgarcx=garcx; tgarcy=garcy! <- Don't switch
944 endif
945 ga=(/tgarcx,tgarcy/)
946 asp=tgarcy/tgarcx
947 basis=basis0
948 
949 call guessak_geo(asp,tgarcx,ak)
950 ma=ga*dtor*0.9_dp ! Shrink first estimate, to start always within bounds
951 
952 ! Perform a Newton iteration (except with imperfect Hessian) to find the
953 ! parameter vector, ak, at which the derivative of Q at constant geographical
954 ! semi-axes, ga, as given, goes to zero. The direct evaluation of the
955 ! Q-derivative at constant ma (which is what is actually computed in
956 ! get_meanq) therefore needs modification to obtain Q-derivarive at constant
957 ! ga:
958 ! dQ/d(ak)|_ga = dQ/d(ak)|_ma - dQ/d(ma)|_ak*d(ma)/d(ga)|_ak*d(ga)/d(ak)|_ma
959 !
960 ! Since the Hessian evaluation is only carried out at constant map-space
961 ! semi-axes, ma, it is not ideal for this problem; consequently, the allowance
962 ! of newton iterations, nit, is much more liberal than we allow for the
963 ! companion routine, bestesg_map, where the constant ma condition WAS
964 ! appropriate.
965 do it=1,nit
966  call get_meanq(ngh,lam,xg,wg,ak,ma,q,qdak,qdma,gat,gadak,gadma,ff)
967  if(ff)return
968  madga=gadma; call inv(madga)! <- d(ma)/d(ga)|_ak ("at constant ak")
969  madak=-matmul(madga,gadak)
970  qdak=qdak+matmul(qdma,madak)! dQ/d(ak)|_ga
971  if(it<=mit)then ! <- Only recompute aks if the basis is new
972 ! Place five additional sample points around the stencil-ellipse:
973  do i=0,4
974  ang=i*dang ! steps of 72 degrees
975  vec2=(/cos(ang),sin(ang)/)*r ! points on a circle of radius r ...
976  dak=matmul(basis,vec2)
977  dma=matmul(madak,dak)
978  aks(:,i)=ak+dak ! ... become points on an ellipse.
979  mas(:,i)=ma+dma
980  enddo
981  call get_meanq(5,ngh,lam,xg,wg,aks,mas, qs,ff)
982  endif
983  grad=matmul(qdak,basis)/q ! <- New grad estimate, accurate to near roundoff
984  if(it<mit)then ! Assume the hessian will have significantly changed
985 ! Recover an approximate Hessian estimate, normalized by q, from all
986 ! 6 samples, q, qs. These are wrt the basis, not (a,k) directly. Note that
987 ! this finite-difference Hessian is wrt q at constant ak, which is roughly
988 ! approximating the desired Hessian wrt q at constant ga; the disparity
989 ! is the reason that we allow more iterations in this routine than in
990 ! companion routine bestesg_map, since we cannot achieve superlinear
991 ! convergence in the present case.
992 ! Make qs the 5-pt discrete Fourier coefficients of the ellipse pts:
993  qs=matmul(em5,qs)/q
994 ! grad=qs(1:2)/r ! Old finite difference estimate is inferior to new grad
995  qs(0)=qs(0)-u1! <- cos(2*ang) coefficient relative to the central value.
996  hess(1,1)=qs(0)+qs(3)! <- combine cos(0*ang) and cos(2*ang) coefficients
997  hess(1,2)=qs(4) ! <- sin(2*ang) coefficient
998  hess(2,1)=qs(4) !
999  hess(2,2)=qs(0)-qs(3)! <- combine cos(0*ang) and cos(2*ang) coefficients
1000  hess=hess*u2/rr ! <- rr is r**2
1001 
1002 ! Perform a Cholesky decomposition of the hessian:
1003  call chol2(hess,el,ff)
1004  if(ff)then
1005  print'(" In bestESG_geo, Hessian is not positive; cholesky fails")'
1006  return
1007  endif
1008 ! Invert the cholesky factor in place:
1009  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)
1010  endif
1011 ! Estimate a Newton step towards the minimum of function Q(a,k):
1012  vec2=-matmul(transpose(el),matmul(el,grad))
1013  dak=matmul(basis,vec2)
1014  gat=gat+matmul(gadak,dak)! <- increment ga
1015  ak=ak+dak! <- increment the parameter vector estimate
1016  dma=-matmul(madga,gat-ga)
1017  ma=ma+dma
1018 
1019 ! Use the inverse cholesky factor to re-condition the basis. This is to make
1020 ! the next stencil-ellipse more closely share the shape of the elliptical
1021 ! contours of Q near its minumum -- essentially a preconditioning of the
1022 ! numerical optimization:
1023  if(it<mit)basis=matmul(basis,transpose(el))
1024 
1025  s=sqrt(dot_product(dak,dak))
1026  if(s<crit)exit ! <-Sufficient convergence of the Newton iteration
1027 enddo
1028 if(it>nit)print'("WARNING; Relatively inferior convergence in bestesg_geo")'
1029 a=ak(1)
1030 k=ak(2)
1031 if(flip)then; marcx=ma(2); marcy=ma(1)! Remember to switch back
1032 else ; marcx=ma(1); marcy=ma(2)! don't switch
1033 endif
1034 end subroutine bestesg_geo
1035 
1072 subroutine bestesg_map(lam,marcx,marcy, a,k,garcx,garcy,q,ff) ![bestesg_map]
1073 use pietc, only: u5,o5,s18,s36,s54,s72,ms18,ms36,ms54,ms72
1074 use psym2, only: chol2
1075 implicit none
1076 real(dp),intent(in ):: lam,marcx,marcy
1077 real(dp),intent(out):: a,k,garcx,garcy,q
1078 logical ,intent(out):: ff
1079 integer(spi),parameter :: nit=25,mit=7,ngh=25
1080 real(dp),parameter :: u2o5=u2*o5, &
1081  f18=u2o5*s18,f36=u2o5*s36,f54=u2o5*s54, &
1082  f72=u2o5*s72,mf18=-f18,mf36=-f36,mf54=-f54,&
1083  mf72=-f72,& !<- (Fourier)
1084  r=0.001_dp,rr=r*r,dang=pi2*o5,crit=1.e-12_dp
1085 real(dp),dimension(ngh) :: wg,xg
1086 real(dp),dimension(0:4,0:4):: em5 ! <- Fourier matrix for 5 points
1087 real(dp),dimension(0:4) :: qs ! <-Sampled q, its Fourier coefficients
1088 real(dp),dimension(2,0:4) :: aks,mas! <- tiny elliptical array of ak
1089 real(dp),dimension(2,2) :: basis0,basis,hess,el,gadak,gadma
1090 real(dp),dimension(2) :: ak,dak,vec2,grad,qdak,qdma,ga,ma
1091 real(dp) :: s,tmarcx,tmarcy,asp,ang
1092 integer(spi) :: i,it
1093 logical :: flip
1094 data em5/o5,u2o5, u0,u2o5, u0,& ! <-The Fourier matrix for 5 points. Applied
1095  o5, f18, f72,mf54, f36,& ! to the five 72-degree spaced values in a
1096  o5,mf54, f36, f18,mf72,& ! column-vector, the product vector has the
1097  o5,mf54,mf36, f18, f72,& ! components, wave-0, cos and sin wave-1,
1098  o5, f18,mf72,mf54,mf36/ ! cos and sin wave-2.
1099 ! First guess upper-triangular basis regularizing the samples used to
1100 ! estimate the Hessian of q by finite differencing:
1101 data basis0/0.1_dp,u0, 0.3_dp,0.3_dp/
1102 ff=lam<u0 .or. lam>=u1
1103 if(ff)then; print'("In bestesg_map; lam out of range")';return; endif
1104 ff= marcx<=u0 .or. marcy<=u0
1105 if(ff)then
1106  print'("In bestesg_map; a nonpositive domain parameter, marcx or marcy")'
1107  return
1108 endif
1109 call gaulegh(ngh,xg,wg)
1110 flip=marcy>marcx
1111 if(flip)then; tmarcx=marcy; tmarcy=marcx! <- Switch
1112 else ; tmarcx=marcx; tmarcy=marcy! <- Don't switch
1113 endif
1114 ma=(/tmarcx,tmarcy/); do i=0,4; mas(:,i)=ma; enddo
1115 asp=tmarcy/tmarcx
1116 basis=basis0
1117 
1118 call guessak_map(asp,tmarcx,ak)
1119 
1120 do it=1,nit
1121  call get_meanq(ngh,lam,xg,wg,ak,ma,q,qdak,qdma,ga,gadak,gadma,ff)
1122  if(ff)return
1123  if(it<=mit)then
1124 ! Place five additional sample points around the stencil-ellipse:
1125  do i=0,4
1126  ang=i*dang ! steps of 72 degrees
1127  vec2=(/cos(ang),sin(ang)/)*r ! points on a circle of radius r ...
1128  aks(:,i)=ak+matmul(basis,vec2) ! ... become points on an ellipse.
1129  enddo
1130  call get_meanq(5,ngh,lam,xg,wg,aks,mas, qs,ff)
1131  endif
1132  grad=matmul(qdak,basis)/q ! <- New grad estimate, accurate to near roundoff
1133  if(it<mit)then
1134 ! Recover Hessian estimate, normalized by q, from all
1135 ! 6 samples, q, qs. These are wrt the basis, not (a,k) directly.
1136 ! The Hessian estimate uses a careful finite method, which is accurate
1137 ! enough. The gradient is NOT estimated by finite differences because we
1138 ! need the gradient to be accurate to near roundoff levels in order that
1139 ! the converged Newton iteration is a precise solution.
1140 ! Make qs the 5-pt discrete Fourier coefficients of the ellipse pts:
1141  qs=matmul(em5,qs)/q
1142 ! grad=qs(1:2)/r ! Old finite difference estimate is inferior to new grad
1143  qs(0)=qs(0)-u1
1144  hess(1,1)=qs(0)+qs(3)! <- combine cos(0*ang) and cos(2*ang) coefficients
1145  hess(1,2)=qs(4) ! <- sin(2*ang) coefficient
1146  hess(2,1)=qs(4) !
1147  hess(2,2)=qs(0)-qs(3)! <- combine cos(0*ang) and cos(2*ang) coefficients
1148  hess=hess*u2/rr ! <- rr is r**2
1149 
1150 ! Perform a Cholesky decomposition of the hessian:
1151  call chol2(hess,el,ff)
1152  if(ff)then
1153  print'(" In bestESG_map, hessian is not positive; cholesky fails")'
1154  return
1155  endif
1156 ! Invert the cholesky factor in place:
1157  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)
1158  endif
1159 ! Estimate a Newton step towards the minimum of function Q(a,k):
1160  vec2=-matmul(transpose(el),matmul(el,grad))
1161  dak=matmul(basis,vec2)
1162  ga=ga+matmul(gadak,dak)! <- increment ga
1163  ak=ak+dak! <- increment the parameter vector estimate
1164 
1165 ! Use the inverse cholesky factor to re-condition the basis. This is to make
1166 ! the next stencil-ellipse more closely share the shape of the elliptical
1167 ! contours of Q near its minumum -- essentially a preconditioning of the
1168 ! numerical optimization:
1169  if(it<mit)basis=matmul(basis,transpose(el))
1170 
1171  s=sqrt(dot_product(dak,dak))
1172  if(s<crit)exit ! <-Sufficient convergence of the Newton iteration
1173 enddo
1174 if(it>nit)print'("WARNING; Relatively inferior convergence in bestesg_map")'
1175 a=ak(1)
1176 k=ak(2)
1177 if(flip)then; garcx=ga(2); garcy=ga(1)! Remember to switch back
1178 else ; garcx=ga(1); garcy=ga(2)! don't switch
1179 endif
1180 end subroutine bestesg_map
1181 
1219 subroutine hgrid_ak_rr(lx,ly,nx,ny,A,K,plat,plon,pazi, & ! [hgrid_ak_rr]
1220  delx,dely, glat,glon,garea, ff)
1221 use pmat4, only: sarea
1222 use pmat5, only: ctogr
1223 implicit none
1224 integer(spi), intent(in ):: lx,ly,nx,ny
1225 real(dp), intent(in ):: a,k,plat,plon,pazi, &
1226  delx,dely
1227 real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: glat,glon
1228 real(dp),dimension(lx:lx+nx-1,ly:ly+ny-1),intent(out):: garea
1229 logical, intent(out):: ff
1230 real(dp),dimension(3,3):: prot,azirot
1231 real(dp),dimension(3,2):: xcd
1232 real(dp),dimension(3) :: xc
1233 real(dp),dimension(2) :: xm
1234 real(dp) :: clat,slat,clon,slon,cazi,sazi,&
1235  rlat,drlata,drlatb,drlatc, &
1236  rlon,drlona,drlonb,drlonc
1237 integer(spi) :: ix,iy,mx,my
1238 clat=cos(plat); slat=sin(plat)
1239 clon=cos(plon); slon=sin(plon)
1240 cazi=cos(pazi); sazi=sin(pazi)
1241 
1242 azirot(:,1)=(/ cazi, sazi, u0/)
1243 azirot(:,2)=(/-sazi, cazi, u0/)
1244 azirot(:,3)=(/ u0, u0, u1/)
1245 
1246 prot(:,1)=(/ -slon, clon, u0/)
1247 prot(:,2)=(/-slat*clon, -slat*slon, clat/)
1248 prot(:,3)=(/ clat*clon, clat*slon, slat/)
1249 prot=matmul(prot,azirot)
1250 mx=lx+nx ! Index of the 'right' edge of the rectangular grid
1251 my=ly+ny ! Index of the 'top' edge of the rectangular grid
1252 do iy=ly,my
1253  xm(2)=iy*dely
1254  do ix=lx,mx
1255  xm(1)=ix*delx
1256  call xmtoxc_ak(a,k,xm,xc,xcd,ff)
1257  if(ff)return
1258  xcd=matmul(prot,xcd)
1259  xc =matmul(prot,xc )
1260  call ctogr(xc,glat(ix,iy),glon(ix,iy))
1261  enddo
1262 enddo
1263 
1264 ! Compute the areas of the quadrilateral grid cells:
1265 do iy=ly,my-1
1266  do ix=lx,mx-1
1267  rlat =glat(ix ,iy )
1268  drlata=glat(ix+1,iy )-rlat
1269  drlatb=glat(ix+1,iy+1)-rlat
1270  drlatc=glat(ix ,iy+1)-rlat
1271  rlon =glon(ix ,iy )
1272  drlona=glon(ix+1,iy )-rlon
1273  drlonb=glon(ix+1,iy+1)-rlon
1274  drlonc=glon(ix ,iy+1)-rlon
1275 ! If 'I' is the grid point (ix,iy), 'A' is (ix+1,iy); 'B' is (ix+1,iy+1)
1276 ! and 'C' is (ix,iy+1), then the area of the grid cell IABC is the sum of
1277 ! the areas of the traingles, IAB and IBC (the latter being the negative
1278 ! of the signed area of triangle, ICB):
1279  garea(ix,iy)=sarea(rlat, drlata,drlona, drlatb,drlonb) &
1280  -sarea(rlat, drlatc,drlonc, drlatb,drlonb)
1281  enddo
1282 enddo
1283 end subroutine hgrid_ak_rr
1284 
1336 subroutine hgrid_ak_rr_c(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak_rr]
1337  delx,dely, glat,glon,garea,dx,dy,angle_dx,angle_dy, ff)
1338 use pmat4, only: cross_product,triple_product
1339 use pmat5, only: ctogr
1340 implicit none
1341 integer(spi), intent(in ):: lx,ly,nx,ny
1342 real(dp), intent(in ):: a,k,plat,plon,pazi, &
1343  delx,dely
1344 real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: glat,glon
1345 real(dp),dimension(lx:lx+nx-1,ly:ly+ny-1),intent(out):: garea
1346 real(dp),dimension(lx:lx+nx-1,ly:ly+ny ),intent(out):: dx
1347 real(dp),dimension(lx:lx+nx ,ly:ly+ny-1),intent(out):: dy
1348 real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: angle_dx,angle_dy
1349 logical, intent(out):: ff
1350 real(dp),dimension(lx-1:lx+nx+1,ly-1:ly+ny+1):: gat ! Temporary area array
1351 real(dp),dimension(lx-1:lx+nx+1,ly :ly+ny ):: dxt ! Temporary dx array
1352 real(dp),dimension(lx :lx+nx ,ly-1:ly+ny+1):: dyt ! Temporary dy array
1353 real(dp),dimension(3,3):: prot,azirot
1354 real(dp),dimension(3,2):: xcd,eano
1355 real(dp),dimension(2,2):: xcd2
1356 real(dp),dimension(3) :: xc,east,north
1357 real(dp),dimension(2) :: xm
1358 real(dp) :: clat,slat,clon,slon,cazi,sazi,delxy
1359 integer(spi) :: ix,iy,mx,my,lxm,lym,mxp,myp
1360 delxy=delx*dely
1361 clat=cos(plat); slat=sin(plat)
1362 clon=cos(plon); slon=sin(plon)
1363 cazi=cos(pazi); sazi=sin(pazi)
1364 
1365 azirot(:,1)=(/ cazi, sazi, u0/)
1366 azirot(:,2)=(/-sazi, cazi, u0/)
1367 azirot(:,3)=(/ u0, u0, u1/)
1368 
1369 prot(:,1)=(/ -slon, clon, u0/)
1370 prot(:,2)=(/-slat*clon, -slat*slon, clat/)
1371 prot(:,3)=(/ clat*clon, clat*slon, slat/)
1372 prot=matmul(prot,azirot)
1373 
1374 mx=lx+nx ! Index of the 'right' edge of the rectangular grid
1375 my=ly+ny ! Index of the 'top' edge of the rectangular grid
1376 lxm=lx-1; mxp=mx+1 ! Indices of extra left and right edges
1377 lym=ly-1; myp=my+1 ! Indices of extra bottom and top edges
1378 
1379 !-- main body of horizontal grid:
1380 do iy=ly,my
1381  xm(2)=iy*dely
1382  do ix=lx,mx
1383  xm(1)=ix*delx
1384  call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
1385  xcd=matmul(prot,xcd)
1386  xc =matmul(prot,xc )
1387  call ctogr(xc,glat(ix,iy),glon(ix,iy))
1388  east=(/-xc(2),xc(1),u0/); east=east/sqrt(dot_product(east,east))
1389  north=cross_product(xc,east)
1390  eano(:,1)=east; eano(:,2)=north
1391  xcd2=matmul(transpose(eano),xcd)
1392  angle_dx(ix,iy)=atan2( xcd2(2,1),xcd2(1,1))
1393  angle_dy(ix,iy)=atan2(-xcd2(1,2),xcd2(2,2))
1394  dxt(ix,iy)=sqrt(dot_product(xcd2(:,1),xcd2(:,1)))*delx
1395  dyt(ix,iy)=sqrt(dot_product(xcd2(:,2),xcd2(:,2)))*dely
1396  gat(ix,iy)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
1397  enddo
1398 enddo
1399 
1400 !-- extra left edge, gat, dxt only:
1401 xm(1)=lxm*delx
1402 do iy=ly,my
1403  xm(2)=iy*dely
1404  call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
1405  xcd=matmul(prot,xcd)
1406  xc =matmul(prot,xc )
1407  east=(/-xc(2),xc(1),u0/); east=east/sqrt(dot_product(east,east))
1408  north=cross_product(xc,east)
1409  eano(:,1)=east; eano(:,2)=north
1410  xcd2=matmul(transpose(eano),xcd)
1411  dxt(lxm,iy)=sqrt(dot_product(xcd2(:,1),xcd2(:,1)))*delx
1412  gat(lxm,iy)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
1413 enddo
1414 
1415 !-- extra right edge, gat, dxt only:
1416 xm(1)=mxp*delx
1417 do iy=ly,my
1418  xm(2)=iy*dely
1419  call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
1420  xcd=matmul(prot,xcd)
1421  xc =matmul(prot,xc )
1422  east=(/-xc(2),xc(1),u0/); east=east/sqrt(dot_product(east,east))
1423  north=cross_product(xc,east)
1424  eano(:,1)=east; eano(:,2)=north
1425  xcd2=matmul(transpose(eano),xcd)
1426  dxt(mxp,iy)=sqrt(dot_product(xcd2(:,1),xcd2(:,1)))*delx
1427  gat(mxp,iy)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
1428 enddo
1429 
1430 !-- extra bottom edge, gat, dyt only:
1431 xm(2)=lym*dely
1432 do ix=lx,mx
1433  xm(1)=ix*delx
1434  call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
1435  xcd=matmul(prot,xcd)
1436  xc =matmul(prot,xc )
1437  east=(/-xc(2),xc(1),u0/); east=east/sqrt(dot_product(east,east))
1438  north=cross_product(xc,east)
1439  eano(:,1)=east; eano(:,2)=north
1440  xcd2=matmul(transpose(eano),xcd)
1441  dyt(ix,lym)=sqrt(dot_product(xcd2(:,2),xcd2(:,2)))*dely
1442  gat(ix,lym)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
1443 enddo
1444 
1445 !-- extra top edge, gat, dyt only:
1446 xm(2)=myp*dely
1447 do ix=lx,mx
1448  xm(1)=ix*delx
1449  call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
1450  xcd=matmul(prot,xcd)
1451  xc =matmul(prot,xc )
1452  east=(/-xc(2),xc(1),u0/); east=east/sqrt(dot_product(east,east))
1453  north=cross_product(xc,east)
1454  eano(:,1)=east; eano(:,2)=north
1455  xcd2=matmul(transpose(eano),xcd)
1456  dyt(ix,myp)=sqrt(dot_product(xcd2(:,2),xcd2(:,2)))*dely
1457  gat(ix,myp)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
1458 enddo
1459 
1460 ! Extra four corners, gat only:
1461 xm(2)=lym*dely
1462 !-- extra bottom left corner:
1463 xm(1)=lxm*delx
1464 call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
1465 xcd=matmul(prot,xcd)
1466 xc =matmul(prot,xc )
1467 gat(lxm,lym)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
1468 
1469 !-- extra bottom right corner:
1470 xm(1)=mxp*delx
1471 call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
1472 xcd=matmul(prot,xcd)
1473 xc =matmul(prot,xc )
1474 gat(mxp,lym)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
1475 
1476 xm(2)=myp*dely
1477 !-- extra top left corner:
1478 xm(1)=lxm*delx
1479 call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
1480 xcd=matmul(prot,xcd)
1481 xc =matmul(prot,xc )
1482 gat(lxm,myp)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
1483 
1484 !-- extra top right corner:
1485 xm(1)=mxp*delx
1486 call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
1487 xcd=matmul(prot,xcd)
1488 xc =matmul(prot,xc )
1489 gat(mxp,myp)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
1490 
1491 !-- 4th-order averaging over each central interval using 4-pt. stencils:
1492 dx =(13_dp*(dxt(lx :mx-1,:)+dxt(lx+1:mx ,:)) &
1493  -(dxt(lxm:mx-2,:)+dxt(lx+2:mxp,:)))/24_dp
1494 dy =(13_dp*(dyt(:,ly :my-1)+dyt(:,ly+1:my )) &
1495  -(dyt(:,lym:my-2)+dyt(:,ly+2:myp)))/24_dp
1496 gat(lx:mx-1,:)=(13_dp*(gat(lx :mx-1,:)+gat(lx+1:mx ,:)) &
1497  -(gat(lxm:mx-2,:)+gat(lx+2:mxp,:)))/24_dp
1498 garea =(13_dp*(gat(lx:mx-1,ly :my-1)+gat(lx:mx-1,ly+1:my )) &
1499  -(gat(lx:mx-1,lym:my-2)+gat(lx:mx-1,ly+2:myp)))/24_dp
1500 end subroutine hgrid_ak_rr_c
1501 
1539 subroutine hgrid_ak_rc(lx,ly,nx,ny,A,K,plat,plon,pazi, & ! [hgrid_ak_rc]
1540  delx,dely, xc,xcd,garea, ff)
1541 use pmat4, only: sarea
1542 use pmat5, only: ctogr
1543 implicit none
1544 integer(spi),intent(in ):: lx,ly,nx,ny
1545 real(dp), intent(in ):: a,k,plat,plon,pazi,delx,dely
1546 real(dp),dimension(3, lx:lx+nx ,ly:ly+ny ),intent(out):: xc
1547 real(dp),dimension(3,2,lx:lx+nx ,ly:ly+ny ),intent(out):: xcd
1548 real(dp),dimension( lx:lx+nx-1,ly:ly+ny-1),intent(out):: garea
1549 logical, intent(out):: ff
1550 real(dp),dimension(3,3):: prot,azirot
1551 real(dp),dimension(2) :: xm
1552 real(dp) :: clat,slat,clon,slon,cazi,sazi, &
1553  rlat,rlata,rlatb,rlatc,drlata,drlatb,drlatc, &
1554  rlon,rlona,rlonb,rlonc,drlona,drlonb,drlonc
1555 integer(spi) :: ix,iy,mx,my
1556 clat=cos(plat); slat=sin(plat)
1557 clon=cos(plon); slon=sin(plon)
1558 cazi=cos(pazi); sazi=sin(pazi)
1559 
1560 azirot(:,1)=(/ cazi, sazi, u0/)
1561 azirot(:,2)=(/-sazi, cazi, u0/)
1562 azirot(:,3)=(/ u0, u0, u1/)
1563 
1564 prot(:,1)=(/ -slon, clon, u0/)
1565 prot(:,2)=(/-slat*clon, -slat*slon, clat/)
1566 prot(:,3)=(/ clat*clon, clat*slon, slat/)
1567 prot=matmul(prot,azirot)
1568 mx=lx+nx ! Index of the 'right' edge of the rectangular grid
1569 my=ly+ny ! Index of the 'top' edge of the rectangular grid
1570 do iy=ly,my
1571  xm(2)=iy*dely
1572  do ix=lx,mx
1573  xm(1)=ix*delx
1574  call xmtoxc_ak(a,k,xm,xc(:,ix,iy),xcd(:,:,ix,iy),ff)
1575  if(ff)return
1576  xcd(:,:,ix,iy)=matmul(prot,xcd(:,:,ix,iy))
1577  xc(:, ix,iy)=matmul(prot,xc(:, ix,iy))
1578  enddo
1579 enddo
1580 
1581 ! Compute the areas of the quadrilateral grid cells:
1582 do iy=ly,my-1
1583  do ix=lx,mx-1
1584  call ctogr(xc(:,ix ,iy ),rlat ,rlon )
1585  call ctogr(xc(:,ix+1,iy ),rlata,rlona)
1586  call ctogr(xc(:,ix+1,iy+1),rlatb,rlonb)
1587  call ctogr(xc(:,ix ,iy+1),rlatc,rlonc)
1588  drlata=rlata-rlat; drlona=rlona-rlon
1589  drlatb=rlatb-rlat; drlonb=rlonb-rlon
1590  drlatc=rlatc-rlat; drlonc=rlonc-rlon
1591 
1592 ! If 'I' is the grid point (ix,iy), 'A' is (ix+1,iy); 'B' is (ix+1,iy+1)
1593 ! and 'C' is (ix,iy+1), then the area of the grid cell IABC is the sum of
1594 ! the areas of the triangles, IAB and IBC (the latter being the negative
1595 ! of the signed area of triangle, ICB):
1596  garea(ix,iy)=sarea(rlat, drlata,drlona, drlatb,drlonb) &
1597  -sarea(rlat, drlatc,drlonc, drlatb,drlonb)
1598  enddo
1599 enddo
1600 end subroutine hgrid_ak_rc
1601 
1628 subroutine hgrid_ak_dd(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, & ! [hgrid_ak_dd]
1629  delx,dely, gdlat,gdlon,garea, ff)
1630 implicit none
1631 integer(spi), intent(in ):: lx,ly,nx,ny
1632 real(dp), intent(in ):: a,k,pdlat,pdlon,&
1633  pdazi,delx,dely
1634 real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: gdlat,gdlon
1635 real(dp),dimension(lx:lx+nx-1,ly:ly+ny-1),intent(out):: garea
1636 logical, intent(out):: ff
1637 real(dp):: plat,plon,pazi
1638 plat=pdlat*dtor ! Convert these angles from degrees to radians
1639 plon=pdlon*dtor !
1640 pazi=pdazi*dtor !
1641 call hgrid_ak_rr(lx,ly,nx,ny,a,k,plat,plon,pazi, &
1642  delx,dely, gdlat,gdlon,garea, ff)
1643 if(ff)return
1644 gdlat=gdlat*rtod ! Convert these angles from radians to degrees
1645 gdlon=gdlon*rtod !
1646 end subroutine hgrid_ak_dd
1647 
1671 subroutine hgrid_ak_dd_c(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, &! [hgrid_ak_dd]
1672  delx,dely, gdlat,gdlon,garea,dx,dy,dangle_dx,dangle_dy, ff)
1673 implicit none
1674 integer(spi), intent(in ):: lx,ly,nx,ny
1675 real(dp), intent(in ):: a,k,pdlat,pdlon,&
1676  pdazi,delx,dely
1677 real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: gdlat,gdlon
1678 real(dp),dimension(lx:lx+nx-1,ly:ly+ny-1),intent(out):: garea
1679 real(dp),dimension(lx:lx+nx-1,ly:ly+ny ),intent(out):: dx
1680 real(dp),dimension(lx:lx+nx ,ly:ly+ny-1),intent(out):: dy
1681 real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: dangle_dx,dangle_dy
1682 logical, intent(out):: ff
1683 real(dp):: plat,plon,pazi
1684 plat=pdlat*dtor ! Convert these angles from degrees to radians
1685 plon=pdlon*dtor !
1686 pazi=pdazi*dtor !
1687 call hgrid_ak_rr_c(lx,ly,nx,ny,a,k,plat,plon,pazi, &
1688  delx,dely, gdlat,gdlon,garea,dx,dy,dangle_dx,dangle_dy, ff)
1689 if(ff)return
1690 gdlat =gdlat *rtod ! Convert these angles from radians to degrees
1691 gdlon =gdlon *rtod !
1692 dangle_dx=dangle_dx*rtod !
1693 dangle_dy=dangle_dy*rtod !
1694 end subroutine hgrid_ak_dd_c
1695 
1722 subroutine hgrid_ak_dc(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, & ! [hgrid_ak_dc]
1723  delx,dely, xc,xcd,garea, ff)
1724 implicit none
1725 integer(spi),intent(in ):: lx,ly,nx,ny
1726 real(dp), intent(in ):: a,k,pdlat,pdlon,pdazi,delx,dely
1727 real(dp),dimension(3, lx:lx+nx ,ly:ly+ny ),intent(out):: xc
1728 real(dp),dimension(3,2,lx:lx+nx ,ly:ly+ny ),intent(out):: xcd
1729 real(dp),dimension( lx:lx+nx-1,ly:ly+ny-1),intent(out):: garea
1730 logical, intent(out):: ff
1731 real(dp):: plat,plon,pazi
1732 plat=pdlat*dtor
1733 plon=pdlon*dtor
1734 pazi=pdazi*dtor
1735 call hgrid_ak_rc(lx,ly,nx,ny,a,k,plat,plon,pazi, &
1736  delx,dely, xc,xcd,garea, ff)
1737 end subroutine hgrid_ak_dc
1738 
1764 subroutine hgrid_ak(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak]
1765  re,delxre,delyre, glat,glon,garea, ff)
1766 implicit none
1767 integer(spi), intent(in ):: lx,ly,nx,ny
1768 real(dp), intent(in ):: a,k,plat,plon,pazi, &
1769  re,delxre,delyre
1770 real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: glat,glon
1771 real(dp),dimension(lx:lx+nx-1,ly:ly+ny-1),intent(out):: garea
1772 logical, intent(out):: ff
1773 real(dp):: delx,dely,rere
1774 delx=delxre/re ! <- set nondimensional grid step delx
1775 dely=delyre/re ! <- set nondimensional grid step dely
1776 call hgrid_ak_rr(lx,ly,nx,ny,a,k,plat,plon,pazi, &
1777  delx,dely, glat,glon,garea, ff)
1778 if(ff)return
1779 rere=re*re
1780 garea=garea*rere ! <- Convert from steradians to physical area units.
1781 end subroutine hgrid_ak
1782 
1816 subroutine hgrid_ak_c(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak]
1817  re,delxre,delyre, glat,glon,garea,dx,dy,dangle_dx,dangle_dy, ff)
1818 implicit none
1819 integer(spi), intent(in ):: lx,ly,nx,ny
1820 real(dp), intent(in ):: a,k,plat,plon,pazi, &
1821  re,delxre,delyre
1822 real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: glat,glon
1823 real(dp),dimension(lx:lx+nx-1,ly:ly+ny-1),intent(out):: garea
1824 real(dp),dimension(lx:lx+nx-1,ly:ly+ny ),intent(out):: dx
1825 real(dp),dimension(lx:lx+nx ,ly:ly+ny-1),intent(out):: dy
1826 real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: dangle_dx,dangle_dy
1827 logical, intent(out):: ff
1828 real(dp):: delx,dely,rere
1829 delx=delxre/re ! <- set nondimensional grid step delx
1830 dely=delyre/re ! <- set nondimensional grid step dely
1831 call hgrid_ak_rr_c(lx,ly,nx,ny,a,k,plat,plon,pazi, &
1832  delx,dely, glat,glon,garea,dx,dy,dangle_dx,dangle_dy, ff)
1833 if(ff)return
1834 rere=re*re
1835 garea=garea*rere ! <- Convert from steradians to physical area units.
1836 dx=dx*re ! <- Convert from nondimensional to physical length units.
1837 dy=dy*re ! <-
1838 dangle_dx=dangle_dx*rtod ! <-Convert from radians to degrees
1839 dangle_dy=dangle_dy*rtod ! <-Convert from radians to degrees
1840 end subroutine hgrid_ak_c
1841 
1851 subroutine gaulegh(m,x,w)! [gaulegh]
1852 implicit none
1853 integer(spi), intent(IN ):: m ! <- number of nodes in half-interval
1854 real(dp),dimension(m),intent(OUT):: x,w ! <- nodes and weights
1855 integer(spi),parameter:: nit=8
1856 real(dp), parameter:: eps=3.e-14_dp
1857 integer(spi) :: i,ic,j,jm,it,m2,m4p,m4p3
1858 real(dp) :: z,zzm,p1,p2,p3,pp,z1
1859 m2=m*2; m4p=m*4+1; m4p3=m4p+2
1860 do i=1,m; ic=m4p3-4*i
1861  z=cos(pih*ic/m4p); zzm=z*z-u1
1862  do it=1,nit
1863  p1=u1; p2=u0
1864  do j=1,m2; jm=j-1; p3=p2; p2=p1; p1=((j+jm)*z*p2-jm*p3)/j; enddo
1865  pp=m2*(z*p1-p2)/zzm; z1=z; z=z1-p1/pp; zzm=z*z-u1
1866  if(abs(z-z1) <= eps)exit
1867  enddo
1868  x(i)=z; w(i)=-u2/(zzm*pp*pp)
1869 enddo
1870 end subroutine gaulegh
1871 
1887 subroutine gtoxm_ak_rr_m(A,K,plat,plon,pazi,lat,lon,xm,ff)! [gtoxm_ak_rr]
1888 use pmat5, only: grtoc
1889 implicit none
1890 real(dp), intent(in ):: a,k,plat,plon,pazi,lat,lon
1891 real(dp),dimension(2),intent(out):: xm
1892 logical, intent(out):: ff
1893 real(dp),dimension(3,3):: prot,azirot
1894 real(dp) :: clat,slat,clon,slon,cazi,sazi
1895 real(dp),dimension(3) :: xc
1896 clat=cos(plat); slat=sin(plat)
1897 clon=cos(plon); slon=sin(plon)
1898 cazi=cos(pazi); sazi=sin(pazi)
1899 
1900 azirot(:,1)=(/ cazi, sazi, u0/)
1901 azirot(:,2)=(/-sazi, cazi, u0/)
1902 azirot(:,3)=(/ u0, u0, u1/)
1903 
1904 prot(:,1)=(/ -slon, clon, u0/)
1905 prot(:,2)=(/-slat*clon, -slat*slon, clat/)
1906 prot(:,3)=(/ clat*clon, clat*slon, slat/)
1907 prot=matmul(prot,azirot)
1908 
1909 call grtoc(lat,lon,xc)
1910 xc=matmul(transpose(prot),xc)
1911 call xctoxm_ak(a,k,xc,xm,ff)
1912 end subroutine gtoxm_ak_rr_m
1913 
1931 subroutine gtoxm_ak_rr_g(A,K,plat,plon,pazi,delx,dely,lat,lon,&! [gtoxm_ak_rr]
1932  xm,ff)
1933 implicit none
1934 real(dp), intent(in ):: a,k,plat,plon,pazi,delx,dely,lat,lon
1935 real(dp),dimension(2),intent(out):: xm
1936 logical, intent(out):: ff
1937 call gtoxm_ak_rr_m(a,k,plat,plon,pazi,lat,lon,xm,ff); if(ff)return
1938 xm(1)=xm(1)/delx; xm(2)=xm(2)/dely
1939 end subroutine gtoxm_ak_rr_g
1940 
1953 subroutine gtoxm_ak_dd_m(A,K,pdlat,pdlon,pdazi,dlat,dlon,&! [gtoxm_ak_dd]
1954  xm,ff)
1955 implicit none
1956 real(dp), intent(in ):: a,k,pdlat,pdlon,pdazi,dlat,dlon
1957 real(dp),dimension(2),intent(out):: xm
1958 logical, intent(out):: ff
1959 real(dp):: plat,plon,pazi,lat,lon
1960 plat=pdlat*dtor ! Convert these angles from degrees to radians
1961 plon=pdlon*dtor !
1962 pazi=pdazi*dtor !
1963 lat=dlat*dtor
1964 lon=dlon*dtor
1965 call gtoxm_ak_rr_m(a,k,plat,plon,pazi,lat,lon,xm,ff)
1966 end subroutine gtoxm_ak_dd_m
1967 
1982 subroutine gtoxm_ak_dd_g(A,K,pdlat,pdlon,pdazi,delx,dely,&! [gtoxm_ak_dd]
1983 dlat,dlon, xm,ff)
1984 implicit none
1985 real(dp), intent(in ):: a,k,pdlat,pdlon,pdazi,delx,dely,dlat,dlon
1986 real(dp),dimension(2),intent(out):: xm
1987 logical, intent(out):: ff
1988 real(dp):: plat,plon,pazi,lat,lon
1989 plat=pdlat*dtor ! Convert these angles from degrees to radians
1990 plon=pdlon*dtor !
1991 pazi=pdazi*dtor !
1992 lat=dlat*dtor
1993 lon=dlon*dtor
1994 call gtoxm_ak_rr_g(a,k,plat,plon,pazi,delx,dely,lat,lon,xm,ff)
1995 end subroutine gtoxm_ak_dd_g
1996 
2014 subroutine xmtog_ak_rr_m(A,K,plat,plon,pazi,xm,lat,lon,ff)! [xmtog_ak_rr]
2015 use pmat5, only: ctogr
2016 implicit none
2017 real(dp), intent(in ):: a,k,plat,plon,pazi
2018 real(dp),dimension(2),intent(in ):: xm
2019 real(dp), intent(out):: lat,lon
2020 logical, intent(out):: ff
2021 real(dp),dimension(3,2):: xcd
2022 real(dp),dimension(3,3):: prot,azirot
2023 real(dp) :: clat,slat,clon,slon,cazi,sazi
2024 real(dp),dimension(3) :: xc
2025 clat=cos(plat); slat=sin(plat)
2026 clon=cos(plon); slon=sin(plon)
2027 cazi=cos(pazi); sazi=sin(pazi)
2028 
2029 azirot(:,1)=(/ cazi, sazi, u0/)
2030 azirot(:,2)=(/-sazi, cazi, u0/)
2031 azirot(:,3)=(/ u0, u0, u1/)
2032 
2033 prot(:,1)=(/ -slon, clon, u0/)
2034 prot(:,2)=(/-slat*clon, -slat*slon, clat/)
2035 prot(:,3)=(/ clat*clon, clat*slon, slat/)
2036 prot=matmul(prot,azirot)
2037 call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
2038 xc=matmul(prot,xc)
2039 call ctogr(xc,lat,lon)
2040 end subroutine xmtog_ak_rr_m
2041 
2061 subroutine xmtog_ak_rr_g(A,K,plat,plon,pazi,delx,dely,xm,&! [xmtog_ak_rr]
2062  lat,lon,ff)
2063 implicit none
2064 real(dp), intent(in ):: a,k,plat,plon,pazi,delx,dely
2065 real(dp),dimension(2),intent(in ):: xm
2066 real(dp), intent(out):: lat,lon
2067 logical, intent(out):: ff
2068 real(dp),dimension(2):: xmt
2069 xmt(1)=xm(1)*delx ! Convert from grid units to intrinsic map-space units
2070 xmt(2)=xm(2)*dely !
2071 call xmtog_ak_rr_m(a,k,plat,plon,pazi,xmt,lat,lon,ff)
2072 end subroutine xmtog_ak_rr_g
2073 
2086 subroutine xmtog_ak_dd_m(A,K,pdlat,pdlon,pdazi,xm,dlat,dlon,ff)! [xmtog_ak_dd]
2087 use pmat5, only: ctogr
2088 implicit none
2089 real(dp), intent(in ):: a,k,pdlat,pdlon,pdazi
2090 real(dp),dimension(2),intent(in ):: xm
2091 real(dp), intent(out):: dlat,dlon
2092 logical, intent(out):: ff
2093 real(dp):: plat,plon,pazi,lat,lon
2094 plat=pdlat*dtor ! Convert these angles from degrees to radians
2095 plon=pdlon*dtor !
2096 pazi=pdazi*dtor !
2097 call xmtog_ak_rr_m(a,k,plat,plon,pazi,xm,lat,lon,ff)
2098 dlat=lat*rtod
2099 dlon=lon*rtod
2100 end subroutine xmtog_ak_dd_m
2101 
2116 subroutine xmtog_ak_dd_g(A,K,pdlat,pdlon,pdazi,delx,dely,xm,&! [xmtog_ak_dd]
2117  dlat,dlon,ff)
2118 implicit none
2119 real(dp), intent(in ):: a,k,pdlat,pdlon,pdazi,delx,dely
2120 real(dp),dimension(2),intent(in ):: xm
2121 real(dp), intent(out):: dlat,dlon
2122 logical, intent(out):: ff
2123 real(dp),dimension(2):: xmt
2124 real(dp) :: plat,plon,pazi,lat,lon
2125 xmt(1)=xm(1)*delx ! Convert from grid units to intrinsic map-space units
2126 xmt(2)=xm(2)*dely !
2127 plat=pdlat*dtor ! Convert these angles from degrees to radians
2128 plon=pdlon*dtor !
2129 pazi=pdazi*dtor !
2130 call xmtog_ak_rr_m(a,k,plat,plon,pazi,xmt,lat,lon,ff)
2131 dlat=lat*rtod
2132 dlon=lon*rtod
2133 end subroutine xmtog_ak_dd_g
2134 
2135 end module pesg
2136 
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:2014
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:602
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:2061
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:770
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:1982
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:1671
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:1887
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:1953
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:2086
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:1816
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:681
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:1931
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:1336
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:2116
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:750
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:541
pesg::gtoxm_ak_rr
Definition: pesg.f90:49