Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_common_coordinate_transformations.F90
Go to the documentation of this file.
1 
106 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
107 #include "sll_assert.h"
108 #include "sll_working_precision.h"
109 
110  use sll_m_constants, only: &
111  sll_p_pi
112 
113  implicit none
114 
115  public :: &
120  sll_f_affine_x1, &
121  sll_f_affine_x2, &
158  sll_f_polar_x1, &
159  sll_f_polar_x2, &
188 
189  private
190 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
191 
192 contains
193 
195  function sll_f_identity_x1(eta1, eta2, params)
196  sll_real64 :: sll_f_identity_x1
197  sll_real64, intent(in) :: eta1
198  sll_real64, intent(in) :: eta2
199  sll_real64, dimension(:), intent(in) :: params
200 #ifdef DEBUG
201  sll_real64 :: dummy
202  dummy = params(1)
203 #endif
204  sll_f_identity_x1 = eta2
205  sll_f_identity_x1 = eta1
206  end function sll_f_identity_x1
207 
209  function sll_f_identity_x2(eta1, eta2, params)
210  sll_real64 :: sll_f_identity_x2
211  sll_real64, intent(in) :: eta1
212  sll_real64, intent(in) :: eta2
213  sll_real64, dimension(:), intent(in) :: params
214 #ifdef DEBUG
215  sll_real64 :: dummy
216  dummy = params(1)
217 #endif
218  sll_f_identity_x2 = eta1
219  sll_f_identity_x2 = eta2
220  end function sll_f_identity_x2
221 
223  function identity_eta1(x1, x2, params)
224  sll_real64 :: identity_eta1
225  sll_real64, intent(in) :: x1
226  sll_real64, intent(in) :: x2
227  sll_real64, dimension(:), intent(in) :: params
228 #ifdef DEBUG
229  sll_real64 :: dummy
230  dummy = params(1)
231 #endif
232  identity_eta1 = x2
233  identity_eta1 = x1
234  end function identity_eta1
235 
237  function identity_eta2(x1, x2, params)
238  sll_real64 :: identity_eta2
239  sll_real64, intent(in) :: x1
240  sll_real64, intent(in) :: x2
241  sll_real64, dimension(:), intent(in) :: params
242 #ifdef DEBUG
243  sll_real64 :: dummy
244  dummy = params(1)
245 #endif
246  identity_eta2 = x1
247  identity_eta2 = x2
248  end function identity_eta2
249 
251  function sll_f_identity_jac11(eta1, eta2, params)
252  sll_real64 :: sll_f_identity_jac11
253  sll_real64, intent(in) :: eta1
254  sll_real64, intent(in) :: eta2
255  sll_real64, dimension(:), intent(in) :: params
256 #ifdef DEBUG
257  sll_real64 :: dummy
258  dummy = params(1)
259 #endif
260  sll_f_identity_jac11 = eta1 + eta2
261  sll_f_identity_jac11 = 1.0_f64
262  end function sll_f_identity_jac11
263 
265  function sll_f_identity_jac12(eta1, eta2, params)
266  sll_real64 :: sll_f_identity_jac12
267  sll_real64, intent(in) :: eta1
268  sll_real64, intent(in) :: eta2
269  sll_real64, dimension(:), intent(in) :: params
270 #ifdef DEBUG
271  sll_real64 :: dummy
272  dummy = params(1)
273 #endif
274  sll_f_identity_jac12 = eta1 + eta2
275  sll_f_identity_jac12 = 0.0_f64
276  end function sll_f_identity_jac12
277 
279  function sll_f_identity_jac21(eta1, eta2, params)
280  sll_real64 :: sll_f_identity_jac21
281  sll_real64, intent(in) :: eta1
282  sll_real64, intent(in) :: eta2
283  sll_real64, dimension(:), intent(in) :: params
284 #ifdef DEBUG
285  sll_real64 :: dummy
286  dummy = params(1)
287 #endif
288  sll_f_identity_jac21 = eta1 + eta2
289  sll_f_identity_jac21 = 0.0_f64
290  end function sll_f_identity_jac21
291 
293  function sll_f_identity_jac22(eta1, eta2, params)
294  sll_real64 :: sll_f_identity_jac22
295  sll_real64, intent(in) :: eta1
296  sll_real64, intent(in) :: eta2
297  sll_real64, dimension(:), intent(in) :: params
298 #ifdef DEBUG
299  sll_real64 :: dummy
300  dummy = params(1)
301 #endif
302  sll_f_identity_jac22 = eta1 + eta2
303  sll_f_identity_jac22 = 1.0_f64
304  end function sll_f_identity_jac22
305 
307  function identity_jac(eta1, eta2, params)
308  sll_real64 :: identity_jac
309  sll_real64, intent(in) :: eta1
310  sll_real64, intent(in) :: eta2
311  sll_real64, dimension(:), intent(in) :: params
312 #ifdef DEBUG
313  sll_real64 :: dummy
314  dummy = params(1)
315 #endif
316  identity_jac = eta1 + eta2
317  identity_jac = 1.0_f64
318  end function identity_jac
319 
321  function sll_f_affine_x1(eta1, eta2, params)
322  sll_real64 :: sll_f_affine_x1
323  sll_real64, intent(in) :: eta1
324  sll_real64, intent(in) :: eta2
325  sll_real64, dimension(:), intent(in) :: params
326  sll_real64 :: a1
327  sll_real64 :: b1
328  sll_real64 :: alpha
329 
330  sll_assert(size(params) >= 5)
331  a1 = params(1)
332  b1 = params(2)
333  alpha = params(5)
334  sll_f_affine_x1 = (b1 - a1)*(cos(alpha)*eta1 - sin(alpha)*eta2) + a1
335  end function sll_f_affine_x1
336 
338  function sll_f_affine_x2(eta1, eta2, params)
339  sll_real64 :: sll_f_affine_x2
340  sll_real64, intent(in) :: eta1
341  sll_real64, intent(in) :: eta2
342  sll_real64, dimension(:), intent(in) :: params
343  sll_real64 :: a2
344  sll_real64 :: b2
345  sll_real64 :: alpha
346 
347  sll_assert(size(params) >= 5)
348  a2 = params(3)
349  b2 = params(4)
350  alpha = params(5)
351  sll_f_affine_x2 = (b2 - a2)*(sin(alpha)*eta1 + cos(alpha)*eta2) + a2
352  end function sll_f_affine_x2
353 
354  ! jacobian maxtrix
356  function sll_f_affine_jac11(eta1, eta2, params)
357  sll_real64 :: sll_f_affine_jac11
358  sll_real64, intent(in) :: eta1
359  sll_real64, intent(in) :: eta2
360  sll_real64, dimension(:), intent(in) :: params
361  sll_real64 :: a1
362  sll_real64 :: b1
363  sll_real64 :: alpha
364 
365  sll_assert(size(params) >= 5)
366  a1 = params(1)
367  b1 = params(2)
368  alpha = params(5)
369  sll_f_affine_jac11 = eta1 + eta2
370  sll_f_affine_jac11 = (b1 - a1)*cos(alpha)
371  end function sll_f_affine_jac11
372 
374  function sll_f_affine_jac12(eta1, eta2, params)
375  sll_real64 :: sll_f_affine_jac12
376  sll_real64, intent(in) :: eta1
377  sll_real64, intent(in) :: eta2
378  sll_real64, dimension(:), intent(in) :: params
379  sll_real64 :: alpha, a1, b1
380  sll_assert(size(params) >= 5)
381  a1 = params(1)
382  b1 = params(2)
383  alpha = params(5)
384  sll_f_affine_jac12 = eta1 + eta2
385  sll_f_affine_jac12 = -(b1 - a1)*sin(alpha)
386  end function sll_f_affine_jac12
387 
389  function sll_f_affine_jac21(eta1, eta2, params)
390  sll_real64 :: sll_f_affine_jac21
391  sll_real64, intent(in) :: eta1
392  sll_real64, intent(in) :: eta2
393  sll_real64 :: alpha, a2, b2
394  sll_real64, dimension(:), intent(in) :: params
395 
396  sll_assert(size(params) >= 5)
397  a2 = params(3)
398  b2 = params(4)
399  alpha = params(5)
400  sll_f_affine_jac21 = eta1 + eta2
401  sll_f_affine_jac21 = (b2 - a2)*sin(alpha)
402  end function sll_f_affine_jac21
403 
405  function sll_f_affine_jac22(eta1, eta2, params)
406  sll_real64 :: sll_f_affine_jac22
407  sll_real64, intent(in) :: eta1
408  sll_real64, intent(in) :: eta2
409  sll_real64, dimension(:), intent(in) :: params
410  sll_real64 :: a2
411  sll_real64 :: b2
412  sll_real64 :: alpha
413 
414  sll_assert(size(params) >= 5)
415  a2 = params(3)
416  b2 = params(4)
417  alpha = params(5)
418  sll_f_affine_jac22 = eta1 + eta2
419  sll_f_affine_jac22 = (b2 - a2)*cos(alpha)
420 
421  end function sll_f_affine_jac22
422 
424  function affine_jac(eta1, eta2, params)
425  sll_real64 :: affine_jac
426  sll_real64, intent(in) :: eta1
427  sll_real64, intent(in) :: eta2
428  sll_real64, dimension(:), intent(in) :: params
429  sll_real64 :: a1
430  sll_real64 :: b1
431  sll_real64 :: a2
432  sll_real64 :: b2
433 
434  sll_assert(size(params) >= 4)
435  a1 = params(1)
436  b1 = params(2)
437  a2 = params(3)
438  b2 = params(4)
439  affine_jac = eta1 + eta2
440  affine_jac = (b1 - a1)*(b2 - a2)
441 
442  end function affine_jac
443 
445  function sll_f_homography_x1(eta1, eta2, params)
446  sll_real64 :: sll_f_homography_x1
447  sll_real64, intent(in) :: eta1
448  sll_real64, intent(in) :: eta2
449  sll_real64, dimension(:), intent(in) :: params
450  sll_real64 :: a, b, c, g, h
451 
452  sll_assert(size(params) >= 8)
453 
454  a = params(1)
455  b = params(2)
456  c = params(3)
457  g = params(7)
458  h = params(8)
459 
460  sll_f_homography_x1 = (a*eta1 + b*eta2 + c)/(g*eta1 + h*eta2 + 1)
461 
462  end function sll_f_homography_x1
463 
465  function sll_f_homography_x2(eta1, eta2, params)
466  sll_real64 :: sll_f_homography_x2
467  sll_real64, intent(in) :: eta1
468  sll_real64, intent(in) :: eta2
469  sll_real64, dimension(:), intent(in) :: params
470  sll_real64 :: d, e, f, g, h
471 
472  sll_assert(size(params) >= 8)
473 
474  d = params(4)
475  e = params(5)
476  f = params(6)
477  g = params(7)
478  h = params(8)
479 
480  sll_f_homography_x2 = (d*eta1 + e*eta2 + f)/(g*eta1 + h*eta2 + 1)
481 
482  end function sll_f_homography_x2
483 
485  function sll_f_homography_jac11(eta1, eta2, params)
486  sll_real64 :: sll_f_homography_jac11
487  sll_real64, intent(in) :: eta1
488  sll_real64, intent(in) :: eta2
489  sll_real64 :: a, b, c, g, h
490  sll_real64, dimension(:), intent(in) :: params
491 
492  sll_assert(size(params) >= 8)
493 
494  a = params(1)
495  b = params(2)
496  c = params(3)
497  g = params(7)
498  h = params(8)
499 
500  sll_f_homography_jac11 = a/(eta1*g + eta2*h + 1) - (a*eta1 + b*eta2 + c)*g/(eta1*g + eta2*h + 1)**2
501 
502  end function sll_f_homography_jac11
503 
505  function sll_f_homography_jac12(eta1, eta2, params)
506  sll_real64 :: sll_f_homography_jac12
507  sll_real64, intent(in) :: eta1
508  sll_real64, intent(in) :: eta2
509  sll_real64, dimension(:), intent(in) :: params
510  sll_real64 :: a, b, c, g, h
511 
512  sll_assert(size(params) >= 8)
513 
514  a = params(1)
515  b = params(2)
516  c = params(3)
517  g = params(7)
518  h = params(8)
519 
520  sll_f_homography_jac12 = b/(eta1*g + eta2*h + 1) - (a*eta1 + b*eta2 + c)*h/(eta1*g + eta2*h + 1)**2
521 
522  end function sll_f_homography_jac12
523 
525  function sll_f_homography_jac21(eta1, eta2, params)
526  sll_real64 :: sll_f_homography_jac21
527  sll_real64, intent(in) :: eta1
528  sll_real64, intent(in) :: eta2
529  sll_real64, dimension(:), intent(in) :: params
530  sll_real64 :: d, e, f, g, h
531 
532  sll_assert(size(params) >= 8)
533 
534  d = params(4)
535  e = params(5)
536  f = params(6)
537  g = params(7)
538  h = params(8)
539 
540  sll_f_homography_jac21 = d/(eta1*g + eta2*h + 1) - (d*eta1 + e*eta2 + f)*g/(eta1*g + eta2*h + 1)**2
541 
542  end function sll_f_homography_jac21
543 
545  function sll_f_homography_jac22(eta1, eta2, params)
546  sll_real64 :: sll_f_homography_jac22
547  sll_real64, intent(in) :: eta1
548  sll_real64, intent(in) :: eta2
549  sll_real64, dimension(:), intent(in) :: params
550  sll_real64 :: d, e, f, g, h
551 
552  sll_assert(size(params) >= 8)
553 
554  d = params(4)
555  e = params(5)
556  f = params(6)
557  g = params(7)
558  h = params(8)
559 
560  sll_f_homography_jac22 = e/(eta1*g + eta2*h + 1) - (d*eta1 + e*eta2 + f)*h/(eta1*g + eta2*h + 1)**2
561 
562  end function sll_f_homography_jac22
563 
565  function homography_jac(eta1, eta2, params)
566  sll_real64 :: homography_jac
567  sll_real64, intent(in) :: eta1
568  sll_real64, intent(in) :: eta2
569  sll_real64, dimension(:), intent(in) :: params
570  sll_real64 :: a, b, c, d, e, f, g, h
571 
572  sll_assert(size(params) >= 8)
573 
574  a = params(1)
575  b = params(2)
576  c = params(3)
577  d = params(4)
578  e = params(5)
579  f = params(6)
580  g = params(7)
581  h = params(8)
582 
583  homography_jac = -(b*d - a*e + (c*e - b*f)*g - (c*d - a*f)*h) &
584  /(eta1**3*g**3 + eta2**3*h**3 + 3*eta1**2*g**2 &
585  + 3*(eta1*eta2**2*g + eta2**2)*h**2 + 3*eta1*g &
586  + 3*(eta1**2*eta2*g**2 + 2*eta1*eta2*g + eta2)*h + 1)
587 
588  end function homography_jac
589 
591  function sll_f_rubber_sheeting_x1(eta1, eta2, params)
592  sll_real64 :: sll_f_rubber_sheeting_x1
593  sll_real64, intent(in) :: eta1
594  sll_real64, intent(in) :: eta2
595  sll_real64, dimension(:), intent(in) :: params
596  sll_real64 :: a, b, c, d
597 
598  sll_assert(size(params) >= 8)
599 
600  a = params(1)
601  b = params(2)
602  c = params(3)
603  d = params(4)
604 
605  sll_f_rubber_sheeting_x1 = a*eta1*eta2 + b*eta1 + c*eta2 + d
606 
607  end function sll_f_rubber_sheeting_x1
608 
610  function sll_f_rubber_sheeting_x2(eta1, eta2, params)
611  sll_real64 :: sll_f_rubber_sheeting_x2
612  sll_real64, intent(in) :: eta1
613  sll_real64, intent(in) :: eta2
614  sll_real64, dimension(:), intent(in) :: params
615  sll_real64 :: e, f, g, h
616 
617  sll_assert(size(params) >= 8)
618 
619  e = params(5)
620  f = params(6)
621  g = params(7)
622  h = params(8)
623 
624  sll_f_rubber_sheeting_x2 = e*eta1*eta2 + f*eta1 + g*eta2 + h
625 
626  end function sll_f_rubber_sheeting_x2
627 
629  function sll_f_rubber_sheeting_jac11(eta1, eta2, params)
630  sll_real64 :: sll_f_rubber_sheeting_jac11
631  sll_real64, intent(in) :: eta1
632  sll_real64, intent(in) :: eta2
633  sll_real64 :: a, b
634  sll_real64, dimension(:), intent(in) :: params
635 
636  sll_assert(size(params) >= 8)
637 
638  a = params(1)
639  b = params(2)
641  sll_f_rubber_sheeting_jac11 = a*eta2 + b
642 
643  end function sll_f_rubber_sheeting_jac11
644 
646  function sll_f_rubber_sheeting_jac12(eta1, eta2, params)
647  sll_real64 :: sll_f_rubber_sheeting_jac12
648  sll_real64, intent(in) :: eta1
649  sll_real64, intent(in) :: eta2
650  sll_real64, dimension(:), intent(in) :: params
651  sll_real64 :: a, c
652 
653  sll_assert(size(params) >= 8)
654 
655  a = params(1)
656  c = params(3)
658  sll_f_rubber_sheeting_jac12 = a*eta1 + c
659 
660  end function sll_f_rubber_sheeting_jac12
661 
663  function sll_f_rubber_sheeting_jac21(eta1, eta2, params)
664  sll_real64 :: sll_f_rubber_sheeting_jac21
665  sll_real64, intent(in) :: eta1
666  sll_real64, intent(in) :: eta2
667  sll_real64, dimension(:), intent(in) :: params
668  sll_real64 :: e, f
669 
670  sll_assert(size(params) >= 8)
671 
672  e = params(5)
673  f = params(6)
675  sll_f_rubber_sheeting_jac21 = e*eta2 + f
676 
677  end function sll_f_rubber_sheeting_jac21
678 
680  function sll_f_rubber_sheeting_jac22(eta1, eta2, params)
681  sll_real64 :: sll_f_rubber_sheeting_jac22
682  sll_real64, intent(in) :: eta1
683  sll_real64, intent(in) :: eta2
684  sll_real64, dimension(:), intent(in) :: params
685  sll_real64 :: e, g
686 
687  sll_assert(size(params) >= 8)
688 
689  e = params(5)
690  g = params(7)
692  sll_f_rubber_sheeting_jac22 = e*eta1 + g
693 
694  end function sll_f_rubber_sheeting_jac22
695 
697  function rubber_sheeting_jac(eta1, eta2, params)
698  sll_real64 :: rubber_sheeting_jac
699  sll_real64, intent(in) :: eta1
700  sll_real64, intent(in) :: eta2
701  sll_real64, dimension(:), intent(in) :: params
702  sll_real64 :: a, b, c, e, f, g
703 
704  sll_assert(size(params) >= 8)
705 
706  a = params(1)
707  b = params(2)
708  c = params(3)
709  e = params(5)
710  f = params(6)
711  g = params(7)
712 
713  rubber_sheeting_jac = (e*eta1 + g)*(a*eta2 + b) - (a*eta1 + c)*(e*eta2 + f)
714 
715  end function rubber_sheeting_jac
716 
717  ! **************************************************************************
718  !
719  ! polar coordinate transformation (r = eta1, theta = eta2):
720  !
721  ! x1 = eta1 * cos (eta2)
722  ! x2 = eta1 * sin (eta2)
723  !
724  ! **************************************************************************
725 
727  function sll_f_polar_x1(eta1, eta2, params)
728  sll_real64 :: sll_f_polar_x1
729  sll_real64, intent(in) :: eta1
730  sll_real64, intent(in) :: eta2
731  sll_real64, dimension(:), intent(in) :: params
732 #ifdef DEBUG
733  sll_real64 :: dummy
734  dummy = params(1)
735 #endif
736  sll_f_polar_x1 = eta1*cos(eta2)
737  end function sll_f_polar_x1
738 
740  function sll_f_polar_x2(eta1, eta2, params)
741  sll_real64 :: sll_f_polar_x2
742  sll_real64, intent(in) :: eta1
743  sll_real64, intent(in) :: eta2
744  sll_real64, dimension(:), intent(in) :: params
745 #ifdef DEBUG
746  sll_real64 :: dummy
747  dummy = params(1)
748 #endif
749  sll_f_polar_x2 = eta1*sin(eta2)
750  end function sll_f_polar_x2
751 
753  function sll_f_polar_eta1(x1, x2, params)
754  sll_real64 :: sll_f_polar_eta1
755  sll_real64, intent(in) :: x1
756  sll_real64, intent(in) :: x2
757  sll_real64, dimension(:), intent(in) :: params
758 #ifdef DEBUG
759  sll_real64 :: dummy
760  dummy = params(1)
761 #endif
762  sll_f_polar_eta1 = sqrt(x1*x1 + x2*x2)
763  end function sll_f_polar_eta1
764 
766  function sll_f_polar_eta2(x1, x2, params)
767  sll_real64 :: sll_f_polar_eta2
768  sll_real64, intent(in) :: x1
769  sll_real64, intent(in) :: x2
770  sll_real64, dimension(:), intent(in) :: params
771 #ifdef DEBUG
772  sll_real64 :: dummy
773  dummy = params(1)
774 #endif
775  sll_f_polar_eta2 = atan(x2/x1)
776  end function sll_f_polar_eta2
777 
779  function sll_f_polar_jac11(eta1, eta2, params)
780  sll_real64 :: sll_f_polar_jac11
781  sll_real64, intent(in) :: eta1
782  sll_real64, intent(in) :: eta2
783  sll_real64, dimension(:), intent(in) :: params
784 #ifdef DEBUG
785  sll_real64 :: dummy
786  dummy = params(1)
787 #endif
788  sll_f_polar_jac11 = eta1
789  sll_f_polar_jac11 = cos(eta2)
790  end function sll_f_polar_jac11
791 
793  function sll_f_polar_jac12(eta1, eta2, params)
794  sll_real64 :: sll_f_polar_jac12
795  sll_real64, intent(in) :: eta1
796  sll_real64, intent(in) :: eta2
797  sll_real64, dimension(:), intent(in) :: params
798 #ifdef DEBUG
799  sll_real64 :: dummy
800  dummy = params(1)
801 #endif
802  sll_f_polar_jac12 = -eta1*sin(eta2)
803  end function sll_f_polar_jac12
804 
806  function sll_f_polar_jac21(eta1, eta2, params)
807  sll_real64 :: sll_f_polar_jac21
808  sll_real64, intent(in) :: eta1
809  sll_real64, intent(in) :: eta2
810  sll_real64, dimension(:), intent(in) :: params
811 #ifdef DEBUG
812  sll_real64 :: dummy
813  dummy = params(1)
814 #endif
815  sll_f_polar_jac21 = eta1
816  sll_f_polar_jac21 = sin(eta2)
817  end function sll_f_polar_jac21
818 
820  function sll_f_polar_jac22(eta1, eta2, params)
821  sll_real64 :: sll_f_polar_jac22
822  sll_real64, intent(in) :: eta1
823  sll_real64, intent(in) :: eta2
824  sll_real64, dimension(:), intent(in) :: params
825 #ifdef DEBUG
826  sll_real64 :: dummy
827  dummy = params(1)
828 #endif
829  sll_f_polar_jac22 = eta1*cos(eta2)
830  end function sll_f_polar_jac22
831 
833  function polar_jac(eta1, eta2, params)
834  sll_real64 :: polar_jac
835  sll_real64, intent(in) :: eta1
836  sll_real64, intent(in) :: eta2
837  sll_real64, dimension(:), intent(in) :: params
838 #ifdef DEBUG
839  sll_real64 :: dummy
840  dummy = params(1)
841 #endif
842  polar_jac = eta2
843  polar_jac = eta1
844  end function polar_jac
845 
846  ! **************************************************************************
847  !
848  ! polar_shear coordinate transformation (r = eta1, theta = eta2):
849  !
850  ! x1 = eta1 * cos (eta2+a*eta1+b)
851  ! x2 = eta1 * sin (eta2+a*eta1+b)
852  !
853  ! **************************************************************************
854 
856  function sll_f_polar_shear_x1(eta1, eta2, params)
857  sll_real64 :: sll_f_polar_shear_x1
858  sll_real64, intent(in) :: eta1
859  sll_real64, intent(in) :: eta2
860  sll_real64, dimension(:), intent(in) :: params
861 
862  sll_real64 :: alpha1
863  sll_real64 :: alpha2
864  alpha1 = params(1)
865  alpha2 = params(2)
866 
867  sll_f_polar_shear_x1 = eta1*cos(eta2 + alpha1*eta1 + alpha2)
868  end function sll_f_polar_shear_x1
869 
871  function sll_f_polar_shear_x2(eta1, eta2, params)
872  sll_real64 :: sll_f_polar_shear_x2
873  sll_real64, intent(in) :: eta1
874  sll_real64, intent(in) :: eta2
875  sll_real64, dimension(:), intent(in) :: params
876 
877  sll_real64 :: alpha1
878  sll_real64 :: alpha2
879  alpha1 = params(1)
880  alpha2 = params(2)
881 
882  sll_f_polar_shear_x2 = eta1*sin(eta2 + alpha1*eta1 + alpha2)
883  end function sll_f_polar_shear_x2
884 
886  function sll_f_polar_shear_jac11(eta1, eta2, params)
887  sll_real64 :: sll_f_polar_shear_jac11
888  sll_real64, intent(in) :: eta1
889  sll_real64, intent(in) :: eta2
890  sll_real64, dimension(:), intent(in) :: params
891 
892  sll_real64 :: alpha1
893  sll_real64 :: alpha2
894  alpha1 = params(1)
895  alpha2 = params(2)
896 
897  sll_f_polar_shear_jac11 = cos(eta2 + alpha1*eta1 + alpha2) &
898  - eta1*alpha1*sin(eta2 + alpha1*eta1 + alpha2)
899  end function sll_f_polar_shear_jac11
900 
902  function sll_f_polar_shear_jac12(eta1, eta2, params)
903  sll_real64 :: sll_f_polar_shear_jac12
904  sll_real64, intent(in) :: eta1
905  sll_real64, intent(in) :: eta2
906  sll_real64, dimension(:), intent(in) :: params
907 
908  sll_real64 :: alpha1
909  sll_real64 :: alpha2
910  alpha1 = params(1)
911  alpha2 = params(2)
912 
913  sll_f_polar_shear_jac12 = -eta1*sin(eta2 + alpha1*eta1 + alpha2)
914  end function sll_f_polar_shear_jac12
915 
917  function sll_f_polar_shear_jac21(eta1, eta2, params)
918  sll_real64 :: sll_f_polar_shear_jac21
919  sll_real64, intent(in) :: eta1
920  sll_real64, intent(in) :: eta2
921  sll_real64, dimension(:), intent(in) :: params
922 
923  sll_real64 :: alpha1
924  sll_real64 :: alpha2
925  alpha1 = params(1)
926  alpha2 = params(2)
927 
928  sll_f_polar_shear_jac21 = sin(eta2 + alpha1*eta1 + alpha2) &
929  + eta1*alpha1*cos(eta2 + alpha1*eta1 + alpha2)
930  end function sll_f_polar_shear_jac21
931 
933  function sll_f_polar_shear_jac22(eta1, eta2, params)
934  sll_real64 :: sll_f_polar_shear_jac22
935  sll_real64, intent(in) :: eta1
936  sll_real64, intent(in) :: eta2
937  sll_real64, dimension(:), intent(in) :: params
938 
939  sll_real64 :: alpha1
940  sll_real64 :: alpha2
941  alpha1 = params(1)
942  alpha2 = params(2)
943 
944  sll_f_polar_shear_jac22 = eta1*cos(eta2 + alpha1*eta1 + alpha2)
945  end function sll_f_polar_shear_jac22
946 
948  function polar_shear_jac(eta1, eta2, params)
949  sll_real64 :: polar_shear_jac
950  sll_real64, intent(in) :: eta1
951  sll_real64, intent(in) :: eta2
952  sll_real64, dimension(:), intent(in) :: params
953 #ifdef DEBUG
954  sll_real64 :: dummy
955  dummy = params(1)
956 #endif
957  polar_shear_jac = eta2
958  polar_shear_jac = eta1
959  end function polar_shear_jac
960 
961  ! **************************************************************************
962  !
963  ! polygonal coordinate transformation (r = eta1, theta = eta2):
964  !
965  ! x1 = eta1 * cos (eta2) * cos(Pi/num_sides) / cos(Pi/num_sides-|eta2|)
966  ! x2 = eta1 * sin (eta2) * cos(Pi/num_sides) / cos(Pi/num_sides-|eta2|)
967  ! for 0 <= eta2 <= 2Pi/num_sides
968  !
969  ! **************************************************************************
970 
972  function sll_f_polygonal_x1(eta1, eta2, params)
973  sll_real64 :: sll_f_polygonal_x1
974  sll_real64, intent(in) :: eta1
975  sll_real64, intent(in) :: eta2
976  sll_real64, dimension(:), intent(in) :: params
977 
978  sll_real64 :: eta2_shift
979  sll_real64 :: eta2r
980  sll_int32 :: eta2i
981  sll_real64 :: eps
982  sll_real64 :: num_sides
983 
984  eps = params(1)
985  num_sides = params(2)
986  eta2_shift = eta2 + params(3)*2._f64*sll_p_pi/num_sides
987  eta2r = eta2 + sll_p_pi/num_sides
988  eta2r = eta2r*0.5_f64*num_sides/sll_p_pi
989  eta2i = floor(eta2r)
990  eta2r = sqrt((eta2 - eta2i*sll_p_pi/(0.5_f64*num_sides))**2 + eps**2)
991  sll_f_polygonal_x1 = eta1*cos(sll_p_pi/num_sides)*cos(eta2_shift)/cos(sll_p_pi/num_sides - eta2r)
992  end function sll_f_polygonal_x1
993 
995  function sll_f_polygonal_x2(eta1, eta2, params)
996  sll_real64 :: sll_f_polygonal_x2
997  sll_real64, intent(in) :: eta1
998  sll_real64, intent(in) :: eta2
999  sll_real64, dimension(:), intent(in) :: params
1000 
1001  sll_real64 :: eta2_shift
1002  sll_real64 :: eta2r
1003  sll_int32 :: eta2i
1004  sll_real64 :: eps
1005  sll_real64 :: num_sides
1006 
1007  eps = params(1)
1008  num_sides = params(2)
1009  eta2_shift = eta2 + params(3)*2._f64*sll_p_pi/num_sides
1010  eta2r = eta2 + sll_p_pi/num_sides
1011  eta2r = eta2r*0.5_f64*num_sides/sll_p_pi
1012  eta2i = floor(eta2r)
1013  eta2r = sqrt((eta2 - eta2i*sll_p_pi/(0.5_f64*num_sides))**2 + eps**2)
1014  sll_f_polygonal_x2 = eta1*cos(sll_p_pi/num_sides)*sin(eta2_shift)/cos(sll_p_pi/num_sides - eta2r)
1015  end function sll_f_polygonal_x2
1016 
1018  function sll_f_polygonal_jac11(eta1, eta2, params)
1019  sll_real64 :: sll_f_polygonal_jac11
1020  sll_real64, intent(in) :: eta1
1021  sll_real64, intent(in) :: eta2
1022  sll_real64, dimension(:), intent(in) :: params
1023 
1024  sll_real64 :: eta2_shift
1025  sll_real64 :: eta2r
1026  sll_int32 :: eta2i
1027  sll_real64 :: eps
1028  sll_real64 :: num_sides
1029 
1030  eps = params(1)
1031  num_sides = params(2)
1032  eta2_shift = eta2 + params(3)*2._f64*sll_p_pi/num_sides
1033  eta2r = eta2 + sll_p_pi/num_sides
1034  eta2r = eta2r*0.5_f64*num_sides/sll_p_pi
1035  eta2i = floor(eta2r)
1036  eta2r = sqrt((eta2 - eta2i*sll_p_pi/(0.5_f64*num_sides))**2 + eps**2)
1037  sll_f_polygonal_jac11 = eta1
1038  sll_f_polygonal_jac11 = cos(sll_p_pi/num_sides)*cos(eta2_shift)/cos(sll_p_pi/num_sides - eta2r)
1039  end function sll_f_polygonal_jac11
1040 
1042  function sll_f_polygonal_jac12(eta1, eta2, params)
1043  sll_real64 :: sll_f_polygonal_jac12
1044  sll_real64, intent(in) :: eta1
1045  sll_real64, intent(in) :: eta2
1046  sll_real64, dimension(:), intent(in) :: params
1047 
1048  sll_real64 :: eta2_shift
1049  sll_real64 :: eta2r
1050  sll_int32 :: eta2i
1051  sll_real64 :: eps
1052  sll_real64 :: num_sides
1053  sll_real64 :: a
1054 
1055  eps = params(1)
1056  num_sides = params(2)
1057  eta2_shift = eta2 + params(3)*2._f64*sll_p_pi/num_sides
1058  eta2r = eta2 + sll_p_pi/num_sides
1059  eta2r = eta2r*0.5_f64*num_sides/sll_p_pi
1060  eta2i = floor(eta2r)
1061  a = eta2i*sll_p_pi/(0.5_f64*num_sides)
1062  eta2r = sqrt((eta2 - a)**2 + eps**2)
1063  sll_f_polygonal_jac12 = -eta1*cos(sll_p_pi/num_sides)*sin(eta2_shift)/cos(sll_p_pi/num_sides - eta2r) &
1064  + eta1*cos(sll_p_pi/num_sides)*cos(eta2_shift) &
1065  *sin(sll_p_pi/num_sides - eta2r)/cos(sll_p_pi/num_sides - eta2r)**2*(eta2 - a)/eta2r
1066  end function sll_f_polygonal_jac12
1067 
1069  function sll_f_polygonal_jac21(eta1, eta2, params)
1070  sll_real64 :: sll_f_polygonal_jac21
1071  sll_real64, intent(in) :: eta1
1072  sll_real64, intent(in) :: eta2
1073  sll_real64, dimension(:), intent(in) :: params
1074 
1075  sll_real64 :: eta2_shift
1076  sll_real64 :: eta2r
1077  sll_int32 :: eta2i
1078  sll_real64 :: eps
1079  sll_real64 :: num_sides
1080 
1081  eps = params(1)
1082  num_sides = params(2)
1083  eta2_shift = eta2 + params(3)*2._f64*sll_p_pi/num_sides
1084  eta2r = eta2 + sll_p_pi/num_sides
1085  eta2r = eta2r*0.5_f64*num_sides/sll_p_pi
1086  eta2i = floor(eta2r)
1087  eta2r = sqrt((eta2 - eta2i*sll_p_pi/(0.5_f64*num_sides))**2 + eps**2)
1088  sll_f_polygonal_jac21 = eta1
1089  sll_f_polygonal_jac21 = cos(sll_p_pi/num_sides)*sin(eta2_shift)/cos(sll_p_pi/num_sides - eta2r)
1090  end function sll_f_polygonal_jac21
1091 
1093  function sll_f_polygonal_jac22(eta1, eta2, params)
1094  sll_real64 :: sll_f_polygonal_jac22
1095  sll_real64, intent(in) :: eta1
1096  sll_real64, intent(in) :: eta2
1097  sll_real64, dimension(:), intent(in) :: params
1098 
1099  sll_real64 :: eta2_shift
1100  sll_real64 :: eta2r
1101  sll_int32 :: eta2i
1102  sll_real64 :: eps
1103  sll_real64 :: num_sides
1104  sll_real64 :: a
1105 
1106  eps = params(1)
1107  num_sides = params(2)
1108  eta2_shift = eta2 + params(3)*2._f64*sll_p_pi/num_sides
1109 
1110  eta2r = eta2 + sll_p_pi/num_sides
1111  eta2r = eta2r*0.5_f64*num_sides/sll_p_pi
1112  eta2i = floor(eta2r)
1113  a = eta2i*sll_p_pi/(0.5_f64*num_sides)
1114  eta2r = sqrt((eta2 - a)**2 + eps**2)
1115  sll_f_polygonal_jac22 = eta1*cos(sll_p_pi/num_sides)*cos(eta2_shift)/cos(sll_p_pi/num_sides - eta2r) &
1116  + eta1*cos(sll_p_pi/num_sides)*sin(eta2_shift) &
1117  *sin(sll_p_pi/num_sides - eta2r)/cos(sll_p_pi/num_sides - eta2r)**2*(eta2 - a)/eta2r
1118  end function sll_f_polygonal_jac22
1119 
1122  function polygonal_jac(eta1, eta2, params)
1123  sll_real64 :: polygonal_jac
1124  sll_real64, intent(in) :: eta1
1125  sll_real64, intent(in) :: eta2
1126  sll_real64, dimension(:), intent(in) :: params
1127 
1128  sll_int32 :: eta2i
1129  sll_real64 :: eta2r
1130 #ifdef DEBUG
1131  sll_real64 :: dummy
1132  dummy = params(1) + eta1
1133 #endif
1134 
1135  eta2r = eta2 + sll_p_pi/6._f64
1136  eta2r = eta2r*3._f64/sll_p_pi
1137  eta2i = floor(eta2r)
1138  eta2r = abs(eta2 - eta2i*sll_p_pi/3._f64)
1139  polygonal_jac = eta2
1140 ! polygonal_jac = sll_f_polygonal_jac11( eta1, eta2, params ) * &
1141 ! sll_f_polygonal_jac22( eta1, eta2, params ) &
1142 ! -sll_f_polygonal_jac21( eta1, eta2, params ) * &
1143 ! sll_f_polygonal_jac12( eta1, eta2, params ) &
1144  end function polygonal_jac
1145 
1146  ! **************************************************************************
1147  !
1148  ! "Colella transformation";
1149  ! sinusoidal product (see P. Colella et al. JCP 230 (2011) formula
1150  ! (102) p 2968):
1151  !
1152  ! x1 = eta1 + alpha1 * sin(2*pi/L1 * eta1) * sin(2*pi/L2 * eta2)
1153  ! x2 = eta2 + alpha2 * sin(2*pi/L1 * eta1) * sin(2*pi/L2 * eta2)
1154  !
1155  ! Domain: [0,L1] X [0,L2]
1156  ! But we generalize it by taking a transformation
1157  ! Domain: [a,b] X [c,d] -----> Domain: [a',b'] X [c',d']
1158  ! so the transformation becomes
1159  !
1160  ! x1 = (b' - a')/(b- a) * ( eta1 + alpha1 * sin(2*pi*( eta1-a)/( b-a)) * sin(2*pi*( eta2-c)/( d-c)) )
1161  ! + ( a' b - b' a)/(b- a)
1162  ! x2 = (d' - c')/(d- c) * ( eta2 + alpha2 * sin(2*pi*( eta1-a)/( b-a)) * sin(2*pi*( eta2-c)/( d-c)) )
1163  ! + ( c' d - d' c)/(d- c)
1164  !
1165  ! The parameters are:
1166  ! alpha1 = params(1)
1167  ! alpha2 = params(2)
1168  ! a = params(3)
1169  ! b = params(4)
1170  ! a' = params(5)
1171  ! b' = params(6)
1172  ! c = params(7)
1173  ! d = params(8)
1174  ! c' = params(9)
1175  ! d' = params(10)
1176  !
1177  !
1178  ! **************************************************************************
1179 
1181  function sll_f_sinprod_gen_x1(eta1, eta2, params)
1182  sll_real64 :: sll_f_sinprod_gen_x1
1183  sll_real64, intent(in) :: eta1
1184  sll_real64, intent(in) :: eta2
1185  sll_real64, dimension(:), intent(in) :: params
1186  sll_real64 :: alpha1
1187  sll_real64 :: alpha2
1188  sll_real64 :: a_1
1189  sll_real64 :: b_1
1190  sll_real64 :: a_2
1191  sll_real64 :: b_2
1192  sll_real64 :: c_1
1193  sll_real64 :: d_1
1194  sll_real64 :: c_2
1195  sll_real64 :: d_2
1196  sll_real64 :: rl1 ! reciprocal of the length of the domain
1197  sll_real64 :: rl2
1198  sll_real64 :: pi2
1199  sll_real64 :: coef1, coef2
1200 
1201  sll_assert(size(params) >= 10)
1202  alpha1 = params(1)
1203  alpha2 = params(2)
1204  a_1 = params(3)
1205  b_1 = params(4)
1206  a_2 = params(5)
1207  b_2 = params(6)
1208  c_1 = params(7)
1209  d_1 = params(8)
1210  c_2 = params(9)
1211  d_2 = params(10)
1212 
1213  rl1 = (eta1 - a_1)/(b_1 - a_1)
1214  rl2 = (eta2 - c_1)/(d_1 - c_1)
1215 
1216  coef1 = (b_2 - a_2)/(b_1 - a_1)
1217  coef2 = (a_2*b_1 - b_2*a_1)/(b_1 - a_1)
1218 
1219  pi2 = 2.0_f64*sll_p_pi
1220  sll_f_sinprod_gen_x1 = eta1 + alpha1*sin(pi2*rl1)*sin(pi2*rl2)
1222  end function sll_f_sinprod_gen_x1
1223 
1225  function sll_f_sinprod_gen_x2(eta1, eta2, params)
1226  sll_real64 :: sll_f_sinprod_gen_x2
1227  sll_real64, intent(in) :: eta1
1228  sll_real64, intent(in) :: eta2
1229  sll_real64, dimension(:), intent(in) :: params
1230  sll_real64 :: alpha1
1231  sll_real64 :: alpha2
1232  sll_real64 :: a_1
1233  sll_real64 :: b_1
1234  sll_real64 :: a_2
1235  sll_real64 :: b_2
1236  sll_real64 :: c_1
1237  sll_real64 :: d_1
1238  sll_real64 :: c_2
1239  sll_real64 :: d_2
1240  sll_real64 :: coef1, coef2
1241  sll_real64 :: rl1 ! reciprocal of the length of the domain
1242  sll_real64 :: rl2
1243  sll_real64 :: pi2
1244 
1245  sll_assert(size(params) >= 10)
1246  alpha1 = params(1)
1247  alpha2 = params(2)
1248  a_1 = params(3)
1249  b_1 = params(4)
1250  a_2 = params(5)
1251  b_2 = params(6)
1252  c_1 = params(7)
1253  d_1 = params(8)
1254  c_2 = params(9)
1255  d_2 = params(10)
1256 
1257  rl1 = (eta1 - a_1)/(b_1 - a_1)
1258  rl2 = (eta2 - c_1)/(d_1 - c_1)
1259 
1260  coef1 = (d_2 - c_2)/(d_1 - c_1)
1261  coef2 = (c_2*d_1 - d_2*c_1)/(d_1 - c_1)
1262 
1263  pi2 = 2.0_f64*sll_p_pi
1264  sll_f_sinprod_gen_x2 = eta2 + alpha2*sin(pi2*rl1)*sin(pi2*rl2)
1266  end function sll_f_sinprod_gen_x2
1267 
1268  ! KK: Unused functions with forcheck errors
1269 !!$ !> inverse mapping
1270 !!$ !> cannot be computed analytically in this case. Use fixed point iterations.
1271 !!$ function sinprod_gen_eta1 ( x1, x2, params )
1272 !!$ sll_real64 :: sinprod_gen_eta1
1273 !!$ sll_real64, intent(in) :: x1
1274 !!$ sll_real64, intent(in) :: x2
1275 !!$ sll_real64, dimension(:), optional, intent(in) :: params
1276 !!$#ifdef DEBUG
1277 !!$ sll_real64 :: dummy
1278 !!$ dummy = params(1)
1279 !!$#endif
1280 !!$ ! NEEDS TO BE IMPLEMENTED
1281 !!$ STOP 'function not implemented'
1282 !!$ sinprod_gen_eta1 = x2
1283 !!$ sinprod_gen_eta1 = x1
1284 !!$ end function sinprod_gen_eta1
1285 !!$
1286 !!$ !> inverse mapping
1287 !!$ !> cannot be computed analytically in this case. Use fixed point iterations.
1288 !!$ function sinprod_gen_eta2 ( x1, x2, params )
1289 !!$ sll_real64 :: sinprod_gen_eta2
1290 !!$ sll_real64, intent(in) :: x1
1291 !!$ sll_real64, intent(in) :: x2
1292 !!$ sll_real64, dimension(:), optional, intent(in) :: params
1293 !!$#ifdef DEBUG
1294 !!$ sll_real64 :: dummy
1295 !!$ dummy = params(1)
1296 !!$#endif
1297 !!$ ! NEEDS TO BE IMPLEMENTED
1298 !!$ STOP 'function not implemented'
1299 !!$ sinprod_gen_eta2 = x1
1300 !!$ sinprod_gen_eta2 = x2
1301 !!$ end function sinprod_gen_eta2
1302 
1304  function sll_f_sinprod_gen_jac11(eta1, eta2, params)
1305  sll_real64 :: sll_f_sinprod_gen_jac11
1306  sll_real64, intent(in) :: eta1
1307  sll_real64, intent(in) :: eta2
1308  sll_real64, dimension(:), intent(in) :: params
1309  sll_real64 :: alpha1
1310  sll_real64 :: alpha2
1311  sll_real64 :: a_1
1312  sll_real64 :: b_1
1313  sll_real64 :: a_2
1314  sll_real64 :: b_2
1315  sll_real64 :: c_1
1316  sll_real64 :: d_1
1317  sll_real64 :: c_2
1318  sll_real64 :: d_2
1319  sll_real64 :: coef1, coef2
1320  sll_real64 :: rl1 ! reciprocal of the length of the domain
1321  sll_real64 :: rl2
1322  sll_real64 :: pi2
1323 
1324  sll_assert(size(params) >= 10)
1325  alpha1 = params(1)
1326  alpha2 = params(2)
1327  a_1 = params(3)
1328  b_1 = params(4)
1329  a_2 = params(5)
1330  b_2 = params(6)
1331  c_1 = params(7)
1332  d_1 = params(8)
1333  c_2 = params(9)
1334  d_2 = params(10)
1335 
1336  rl1 = (eta1 - a_1)/(b_1 - a_1)
1337  rl2 = (eta2 - c_1)/(d_1 - c_1)
1338 
1339  coef1 = (b_2 - a_2)/(b_1 - a_1)
1340  coef2 = (a_2*b_1 - b_2*a_1)/(b_1 - a_1)
1341 
1342  pi2 = 2.0_f64*sll_p_pi
1343 
1344  sll_f_sinprod_gen_jac11 = 1.0_f64 + alpha1*pi2/(b_1 - a_1)*cos(pi2*rl1)*sin(pi2*rl2)
1346  end function sll_f_sinprod_gen_jac11
1347 
1349  function sll_f_sinprod_gen_jac12(eta1, eta2, params)
1350  sll_real64 :: sll_f_sinprod_gen_jac12
1351  sll_real64, intent(in) :: eta1
1352  sll_real64, intent(in) :: eta2
1353  sll_real64, dimension(:), intent(in) :: params
1354  sll_real64 :: alpha1
1355  sll_real64 :: alpha2
1356  sll_real64 :: a_1
1357  sll_real64 :: b_1
1358  sll_real64 :: a_2
1359  sll_real64 :: b_2
1360  sll_real64 :: c_1
1361  sll_real64 :: d_1
1362  sll_real64 :: c_2
1363  sll_real64 :: d_2
1364  sll_real64 :: coef1, coef2
1365  sll_real64 :: rl1 ! reciprocal of the length of the domain
1366  sll_real64 :: rl2
1367  sll_real64 :: pi2
1368 
1369  sll_assert(size(params) >= 10)
1370  alpha1 = params(1)
1371  alpha2 = params(2)
1372  a_1 = params(3)
1373  b_1 = params(4)
1374  a_2 = params(5)
1375  b_2 = params(6)
1376  c_1 = params(7)
1377  d_1 = params(8)
1378  c_2 = params(9)
1379  d_2 = params(10)
1380 
1381  rl1 = (eta1 - a_1)/(b_1 - a_1)
1382  rl2 = (eta2 - c_1)/(d_1 - c_1)
1383 
1384  coef1 = (b_2 - a_2)/(b_1 - a_1)
1385  coef2 = (a_2*b_1 - b_2*a_1)/(b_1 - a_1)
1386 
1387  pi2 = 2.0_f64*sll_p_pi
1388 
1389  sll_f_sinprod_gen_jac12 = alpha1*pi2/(d_1 - c_1)*sin(pi2*rl1)*cos(pi2*rl2)
1391  end function sll_f_sinprod_gen_jac12
1392 
1394  function sll_f_sinprod_gen_jac21(eta1, eta2, params)
1395  sll_real64 :: sll_f_sinprod_gen_jac21
1396  sll_real64, intent(in) :: eta1
1397  sll_real64, intent(in) :: eta2
1398  sll_real64, dimension(:), intent(in) :: params
1399  sll_real64 :: alpha1
1400  sll_real64 :: alpha2
1401  sll_real64 :: a_1
1402  sll_real64 :: b_1
1403  sll_real64 :: a_2
1404  sll_real64 :: b_2
1405  sll_real64 :: c_1
1406  sll_real64 :: d_1
1407  sll_real64 :: c_2
1408  sll_real64 :: d_2
1409  sll_real64 :: coef1, coef2
1410  sll_real64 :: rl1 ! reciprocal of the length of the domain
1411  sll_real64 :: rl2
1412  sll_real64 :: pi2
1413 
1414  sll_assert(size(params) >= 10)
1415  alpha1 = params(1)
1416  alpha2 = params(2)
1417  a_1 = params(3)
1418  b_1 = params(4)
1419  a_2 = params(5)
1420  b_2 = params(6)
1421  c_1 = params(7)
1422  d_1 = params(8)
1423  c_2 = params(9)
1424  d_2 = params(10)
1425 
1426  rl1 = (eta1 - a_1)/(b_1 - a_1)
1427  rl2 = (eta2 - c_1)/(d_1 - c_1)
1428 
1429  coef1 = (d_2 - c_2)/(d_1 - c_1)
1430  coef2 = (c_2*d_1 - d_2*c_1)/(d_1 - c_1)
1431 
1432  pi2 = 2.0_f64*sll_p_pi
1433  sll_f_sinprod_gen_jac21 = alpha2*pi2/(b_1 - a_1)*cos(pi2*rl1)*sin(pi2*rl2)
1435  end function sll_f_sinprod_gen_jac21
1436 
1438  function sll_f_sinprod_gen_jac22(eta1, eta2, params)
1439  sll_real64 :: sll_f_sinprod_gen_jac22
1440  sll_real64, intent(in) :: eta1
1441  sll_real64, intent(in) :: eta2
1442  sll_real64, dimension(:), intent(in) :: params
1443  sll_real64 :: alpha1
1444  sll_real64 :: alpha2
1445  sll_real64 :: a_1
1446  sll_real64 :: b_1
1447  sll_real64 :: a_2
1448  sll_real64 :: b_2
1449  sll_real64 :: c_1
1450  sll_real64 :: d_1
1451  sll_real64 :: c_2
1452  sll_real64 :: d_2
1453  sll_real64 :: coef1, coef2
1454  sll_real64 :: rl1 ! reciprocal of the length of the domain
1455  sll_real64 :: rl2
1456  sll_real64 :: pi2
1457 
1458  sll_assert(size(params) >= 10)
1459  alpha1 = params(1)
1460  alpha2 = params(2)
1461  a_1 = params(3)
1462  b_1 = params(4)
1463  a_2 = params(5)
1464  b_2 = params(6)
1465  c_1 = params(7)
1466  d_1 = params(8)
1467  c_2 = params(9)
1468  d_2 = params(10)
1469 
1470  rl1 = (eta1 - a_1)/(b_1 - a_1)
1471  rl2 = (eta2 - c_1)/(d_1 - c_1)
1472 
1473  coef1 = (d_2 - c_2)/(d_1 - c_1)
1474  coef2 = (c_2*d_1 - d_2*c_1)/(d_1 - c_1)
1475 
1476  pi2 = 2.0_f64*sll_p_pi
1477 
1478  sll_f_sinprod_gen_jac22 = 1.0_f64 + alpha2*pi2/(d_1 - c_1)*sin(pi2*rl1)*cos(pi2*rl2)
1480 
1481  end function sll_f_sinprod_gen_jac22
1482 
1484  function sll_f_sinprod_gen_jac(eta1, eta2, params)
1485  sll_real64 :: sll_f_sinprod_gen_jac
1486  sll_real64, intent(in) :: eta1
1487  sll_real64, intent(in) :: eta2
1488  sll_real64, dimension(:), intent(in) :: params
1489 
1490  sll_assert(size(params) >= 10)
1491  !alpha1 = params(1)
1492  !alpha2 = params(2)
1493  !rl1 = 1.0_f64/params(3)
1494  !rl2 = 1.0_f64/params(4)
1495  !pi2 = 2.0_f64*sll_p_pi
1496  !sll_f_sinprod_jac = 1.0_f64 + 0.2_f64 *sll_p_pi * sin (2*sll_p_pi**(eta1+eta2))
1497  !sll_f_sinprod_gen_jac = 1.0_f64 + alpha2*pi2*rl2*sin(pi2*rl1*eta1)*cos(pi2*rl2*eta2) + &
1498  ! alpha1*pi2*rl1*cos(pi2*rl1*eta1)*sin(pi2*rl2*eta2)
1499 
1500  sll_f_sinprod_gen_jac = sll_f_sinprod_gen_jac22(eta1, eta2, params)*sll_f_sinprod_gen_jac11(eta1, eta2, params) &
1501  - sll_f_sinprod_gen_jac12(eta1, eta2, params)*sll_f_sinprod_gen_jac21(eta1, eta2, params)
1502  end function sll_f_sinprod_gen_jac
1503 
1504  ! **************************************************************************
1505  !
1506  ! "Colella transformation";
1507  ! sinusoidal product (see P. Colella et al. JCP 230 (2011) formula
1508  ! (102) p 2968):
1509  !
1510  ! x1 = eta1 + alpha1 * sin(2*pi*eta1) * sin(2*pi*eta2)
1511  ! x2 = eta2 + alpha2 * sin(2*pi*eta1) * sin(2*pi*eta2)
1512  !
1513  ! Domain: [0,L1] X [0,L2]
1514  ! By default the values of the alpha parameters are:
1515  ! alpha1 = 0.1
1516  ! alpha2 = 0.1
1517  ! L1 = 1.0
1518  ! L2 = 1.0
1519  !
1520  ! These parameters are stored in the params array as:
1521  ! ( alpha1, alpha2, L1, L2 )
1522  !
1523  ! **************************************************************************
1524 
1526  function sll_f_sinprod_x1(eta1, eta2, params)
1527  sll_real64 :: sll_f_sinprod_x1
1528  sll_real64, intent(in) :: eta1
1529  sll_real64, intent(in) :: eta2
1530  sll_real64, dimension(:), intent(in) :: params
1531  sll_real64 :: alpha1
1532  sll_real64 :: rl1 ! reciprocal of the length of the domain
1533  sll_real64 :: rl2
1534  sll_real64 :: pi2
1535 
1536  sll_assert(size(params) >= 4)
1537  alpha1 = params(1)
1538  rl1 = 1.0_f64/params(3)
1539  rl2 = 1.0_f64/params(4)
1540  pi2 = 2.0_f64*sll_p_pi
1541  sll_f_sinprod_x1 = eta1 + alpha1*sin(pi2*rl1*eta1)*sin(pi2*rl2*eta2)
1542  end function sll_f_sinprod_x1
1543 
1545  function sll_f_sinprod_x2(eta1, eta2, params)
1546  sll_real64 :: sll_f_sinprod_x2
1547  sll_real64, intent(in) :: eta1
1548  sll_real64, intent(in) :: eta2
1549  sll_real64, dimension(:), intent(in) :: params
1550  sll_real64 :: alpha2
1551  sll_real64 :: rl1 ! reciprocal of the length of the domain
1552  sll_real64 :: rl2
1553  sll_real64 :: pi2
1554 
1555  sll_assert(size(params) >= 4)
1556  alpha2 = params(2)
1557  rl1 = 1.0_f64/params(3)
1558  rl2 = 1.0_f64/params(4)
1559  pi2 = 2.0_f64*sll_p_pi
1560  sll_f_sinprod_x2 = eta2 + alpha2*sin(pi2*rl1*eta1)*sin(pi2*rl2*eta2)
1561  end function sll_f_sinprod_x2
1562 
1563 ! KK: Unused functions with forcheck errors
1564 !!$ !> inverse mapping
1565 !!$ !> cannot be computed analytically in this case. Use fixed point iterations.
1566 !!$ function sinprod_eta1 ( x1, x2, params )
1567 !!$ sll_real64 :: sinprod_eta1
1568 !!$ sll_real64, intent(in) :: x1
1569 !!$ sll_real64, intent(in) :: x2
1570 !!$ sll_real64, dimension(:), optional, intent(in) :: params
1571 !!$#ifdef DEBUG
1572 !!$ sll_real64 :: dummy
1573 !!$ dummy = params(1)
1574 !!$#endif
1575 !!$ ! NEEDS TO BE IMPLEMENTED
1576 !!$ STOP 'function not implemented'
1577 !!$ sinprod_eta1 = x2
1578 !!$ sinprod_eta1 = x1
1579 !!$ end function sinprod_eta1
1580 !!$
1581 !!$ !> inverse mapping
1582 !!$ !> cannot be computed analytically in this case. Use fixed point iterations.
1583 !!$ function sinprod_eta2 ( x1, x2, params )
1584 !!$ sll_real64 :: sinprod_eta2
1585 !!$ sll_real64, intent(in) :: x1
1586 !!$ sll_real64, intent(in) :: x2
1587 !!$ sll_real64, dimension(:), optional, intent(in) :: params
1588 !!$#ifdef DEBUG
1589 !!$ sll_real64 :: dummy
1590 !!$ dummy = params(1)
1591 !!$#endif
1592 !!$ ! NEEDS TO BE IMPLEMENTED
1593 !!$ STOP 'function not implemented'
1594 !!$ sinprod_eta2 = x1
1595 !!$ sinprod_eta2 = x2
1596 !!$ end function sinprod_eta2
1597 
1599  function sll_f_sinprod_jac11(eta1, eta2, params)
1600  sll_real64 :: sll_f_sinprod_jac11
1601  sll_real64, intent(in) :: eta1
1602  sll_real64, intent(in) :: eta2
1603  sll_real64, dimension(:), intent(in) :: params
1604  sll_real64 :: alpha1
1605  sll_real64 :: rl1 ! reciprocal of the length of the domain
1606  sll_real64 :: rl2
1607  sll_real64 :: pi2
1608 
1609  sll_assert(size(params) >= 4)
1610  alpha1 = params(1)
1611  rl1 = 1.0_f64/params(3)
1612  rl2 = 1.0_f64/params(4)
1613  pi2 = 2.0_f64*sll_p_pi
1614  sll_f_sinprod_jac11 = 1.0_f64 + alpha1*pi2*rl1*cos(pi2*rl1*eta1)*sin(pi2*rl2*eta2)
1615  end function sll_f_sinprod_jac11
1616 
1618  function sll_f_sinprod_jac12(eta1, eta2, params)
1619  sll_real64 :: sll_f_sinprod_jac12
1620  sll_real64, intent(in) :: eta1
1621  sll_real64, intent(in) :: eta2
1622  sll_real64, dimension(:), intent(in) :: params
1623  sll_real64 :: alpha1
1624  sll_real64 :: rl1 ! reciprocal of the length of the domain
1625  sll_real64 :: rl2
1626  sll_real64 :: pi2
1627 
1628  sll_assert(size(params) >= 4)
1629  alpha1 = params(1)
1630  rl1 = 1.0_f64/params(3)
1631  rl2 = 1.0_f64/params(4)
1632  pi2 = 2.0_f64*sll_p_pi
1633  sll_f_sinprod_jac12 = alpha1*pi2*rl2*sin(pi2*rl1*eta1)*cos(pi2*rl2*eta2)
1634  end function sll_f_sinprod_jac12
1635 
1637  function sll_f_sinprod_jac21(eta1, eta2, params)
1638  sll_real64 :: sll_f_sinprod_jac21
1639  sll_real64, intent(in) :: eta1
1640  sll_real64, intent(in) :: eta2
1641  sll_real64, dimension(:), intent(in) :: params
1642  sll_real64 :: alpha2
1643  sll_real64 :: rl1 ! reciprocal of the length of the domain
1644  sll_real64 :: rl2
1645  sll_real64 :: pi2
1646 
1647  sll_assert(size(params) >= 4)
1648  alpha2 = params(2)
1649  rl1 = 1.0_f64/params(3)
1650  rl2 = 1.0_f64/params(4)
1651  pi2 = 2.0_f64*sll_p_pi
1652  sll_f_sinprod_jac21 = alpha2*pi2*rl1*cos(pi2*rl1*eta1)*sin(pi2*rl2*eta2)
1653  end function sll_f_sinprod_jac21
1654 
1656  function sll_f_sinprod_jac22(eta1, eta2, params)
1657  sll_real64 :: sll_f_sinprod_jac22
1658  sll_real64, intent(in) :: eta1
1659  sll_real64, intent(in) :: eta2
1660  sll_real64, dimension(:), intent(in) :: params
1661  sll_real64 :: alpha2
1662  sll_real64 :: rl1 ! reciprocal of the length of the domain
1663  sll_real64 :: rl2
1664  sll_real64 :: pi2
1665 
1666  sll_assert(size(params) >= 4)
1667  alpha2 = params(2)
1668  rl1 = 1.0_f64/params(3)
1669  rl2 = 1.0_f64/params(4)
1670  pi2 = 2.0_f64*sll_p_pi
1671  sll_f_sinprod_jac22 = 1.0_f64 + alpha2*pi2*rl2*sin(pi2*rl1*eta1)*cos(pi2*rl2*eta2)
1672  end function sll_f_sinprod_jac22
1673 
1675  function sll_f_sinprod_jac(eta1, eta2, params)
1676  sll_real64 :: sll_f_sinprod_jac
1677  sll_real64, intent(in) :: eta1
1678  sll_real64, intent(in) :: eta2
1679  sll_real64, dimension(:), intent(in) :: params
1680  sll_real64 :: alpha1
1681  sll_real64 :: alpha2
1682  sll_real64 :: rl1 ! reciprocal of the length of the domain
1683  sll_real64 :: rl2
1684  sll_real64 :: pi2
1685 
1686  sll_assert(size(params) >= 4)
1687  alpha1 = params(1)
1688  alpha2 = params(2)
1689  rl1 = 1.0_f64/params(3)
1690  rl2 = 1.0_f64/params(4)
1691  pi2 = 2.0_f64*sll_p_pi
1692  !sll_f_sinprod_jac = 1.0_f64 + 0.2_f64 *sll_p_pi * sin (2*sll_p_pi**(eta1+eta2))
1693  sll_f_sinprod_jac = 1.0_f64 + alpha2*pi2*rl2*sin(pi2*rl1*eta1)*cos(pi2*rl2*eta2) + &
1694  alpha1*pi2*rl1*cos(pi2*rl1*eta1)*sin(pi2*rl2*eta2)
1695  end function sll_f_sinprod_jac
1696 
1697 #if 0
1698 ! Only one Colella transformation should survive, the parametrized one above...
1699 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1700 
1701  ! direct mapping collela on (0,4*pi)x (0,1)
1702  ! Same story as above, we separate the alpha coefficients in each direction
1703  ! and give default values of 0.1. These and the previous transformation
1704  ! should be merged and parametrized with params.
1705  function sinprod_x1_rect(eta1, eta2, params)
1706  real(8) :: sinprod_x1_rect
1707  real(8), intent(in) :: eta1
1708  real(8), intent(in) :: eta2
1709  sll_real64, dimension(:), optional, intent(in) :: params
1710  sll_real64 :: alpha1
1711  sll_real64 :: pi2
1712  if (present(params)) then
1713  sll_assert(size(params) >= 2)
1714  alpha1 = params(1)
1715  else
1716  alpha1 = 0.1_f64
1717  end if
1718  pi2 = 2.0_f64*sll_p_pi
1719  sinprod_x1_rect = eta1 + alpha1*sin(0.5*eta1)*sin(pi2*eta2)
1720  end function sinprod_x1_rect
1721 
1722  function sinprod_x2_rect(eta1, eta2, params)
1723  real(8) :: sinprod_x2_rect
1724  real(8), intent(in) :: eta1
1725  real(8), intent(in) :: eta2
1726  sll_real64, dimension(:), optional, intent(in) :: params
1727  sll_real64 :: alpha2
1728  sll_real64 :: pi2
1729  if (present(params)) then
1730  sll_assert(size(params) >= 2)
1731  alpha2 = params(2)
1732  else
1733  alpha2 = 0.1_f64
1734  end if
1735  pi2 = 2.0_f64*sll_p_pi
1736  sinprod_x2_rect = eta2 + alpha2*sin(0.5*eta1)*sin(pi2*eta2)
1737  end function sinprod_x2_rect
1738 
1739  ! jacobian matrix
1740  function sinprod_jac11_rect(eta1, eta2, params)
1741  real(8) :: sinprod_jac11_rect
1742  real(8), intent(in) :: eta1
1743  real(8), intent(in) :: eta2
1744  sll_real64, dimension(:), optional, intent(in) :: params
1745  if (present(params)) then
1746  sll_assert(size(params) >= 2)
1747  alpha1 = params(1)
1748  else
1749  alpha1 = 0.1_f64
1750  end if
1751  sinprod_jac11_rect = 1.0_8 + 0.5_8*alpha1*cos(0.5*eta1)*sin(2*sll_p_pi*eta2)
1752  end function sinprod_jac11_rect
1753 
1754  function sinprod_jac12_rect(eta1, eta2, params)
1755  real(8) :: sinprod_jac12_rect
1756  real(8), intent(in) :: eta1
1757  real(8), intent(in) :: eta2
1758  sll_real64, dimension(:), optional, intent(in) :: params
1759  if (present(params)) then
1760  sll_assert(size(params) >= 2)
1761  alpha1 = params(1)
1762  else
1763  alpha1 = 0.1_f64
1764  end if
1765  sinprod_jac12_rect = 2._8*sll_p_pi*alpha1*sin(0.5*eta1)*cos(2*sll_p_pi*eta2)
1766  end function sinprod_jac12_rect
1767 
1768  function sinprod_jac21_rect(eta1, eta2, params)
1769  real(8) :: sinprod_jac21_rect
1770  real(8), intent(in) :: eta1
1771  real(8), intent(in) :: eta2
1772  sll_real64, dimension(:), optional, intent(in) :: params
1773  if (present(params)) then
1774  sll_assert(size(params) >= 2)
1775  alpha2 = params(2)
1776  else
1777  alpha2 = 0.1_f64
1778  end if
1779  sinprod_jac21_rect = 0.5_8*alpha2*cos(0.5*eta1)*sin(2*sll_p_pi*eta2)
1780  end function sinprod_jac21_rect
1781 
1782  function sinprod_jac22_rect(eta1, eta2, params)
1783  real(8) :: sinprod_jac22_rect
1784  real(8), intent(in) :: eta1
1785  real(8), intent(in) :: eta2
1786  sll_real64, dimension(:), optional, intent(in) :: params
1787  if (present(params)) then
1788  sll_assert(size(params) >= 2)
1789  alpha2 = params(2)
1790  else
1791  alpha2 = 0.1_f64
1792  end if
1793  sinprod_jac22_rect = 1.0_8 + &
1794  2.0_8*sll_p_pi*alpha2*sin(0.5*eta1)*cos(2*sll_p_pi*eta2)
1795  end function sinprod_jac22_rect
1796 
1797  ! jacobian ie determinant of jacobian matrix
1798  function sinprod_jac_rect(eta1, eta2, params)
1799  real(8) :: sinprod_jac_rect
1800  real(8), intent(in) :: eta1
1801  real(8), intent(in) :: eta2
1802  sll_real64, dimension(:), optional, intent(in) :: params
1803  !sll_f_sinprod_jac = 1.0_f64 + 0.2_f64 *sll_p_pi * sin (2*sll_p_pi**(eta1+eta2))
1804  if (present(params)) then
1805  sll_assert(size(params) >= 2)
1806  alpha2 = params(2)
1807  alpha1 = params(1)
1808  else
1809  alpha2 = 0.1_f64
1810  alpha1 = 0.1_f64
1811  end if
1812  sinprod_jac_rect = &
1813  (1.0_8 + 0.5_8*alpha1*cos(0.5*eta1)*sin(2*sll_p_pi*eta2))* &
1814  (1.0_8 + 2.0_8*sll_p_pi*alpha2*sin(0.5*eta1)*cos(2*sll_p_pi*eta2)) - &
1815  2.0_8*sll_p_pi*alpha2*sin(0.5*eta1)*cos(2*sll_p_pi*eta2)* &
1816  0.5_8*alpha1*cos(0.5*eta1)*sin(2*sll_p_pi*eta2)
1817  end function sinprod_jac_rect
1818 
1819 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1820 
1821  ! direct mapping collela on (0,4*pi)x (0,4*pi)
1822  function sinprod_x1_square(eta1, eta2, params)
1823  real(8) :: sinprod_x1_square
1824  real(8), intent(in) :: eta1
1825  real(8), intent(in) :: eta2
1826  sll_real64, dimension(:), optional, intent(in) :: params
1827  sinprod_x1_square = eta1 + 0.1_8*sin(0.5*eta1)*sin(0.5*eta2)
1828  end function sinprod_x1_square
1829 
1830  function sinprod_x2_square(eta1, eta2, params)
1831  real(8) :: sinprod_x2_square
1832  real(8), intent(in) :: eta1
1833  real(8), intent(in) :: eta2
1834  sll_real64, dimension(:), optional, intent(in) :: params
1835  sinprod_x2_square = eta2 + 0.1_8*sin(0.5*eta1)*sin(0.5*eta2)
1836  end function sinprod_x2_square
1837 
1838  ! jacobian matrix
1839  function sinprod_jac11_square(eta1, eta2, params)
1840  real(8) :: sinprod_jac11_square
1841  real(8), intent(in) :: eta1
1842  real(8), intent(in) :: eta2
1843  sll_real64, dimension(:), optional, intent(in) :: params
1844  sinprod_jac11_square = 1.0_8 + &
1845  0.5_8*0.1_8*cos(0.5*eta1)*sin(0.5*eta2)
1846  end function sinprod_jac11_square
1847 
1848  function sinprod_jac12_square(eta1, eta2, params)
1849  real(8) :: sinprod_jac12_square
1850  real(8), intent(in) :: eta1
1851  real(8), intent(in) :: eta2
1852  sll_real64, dimension(:), optional, intent(in) :: params
1853  sinprod_jac12_square = 0.5_8*0.1_8*sin(0.5*eta1)*cos(0.5*eta2)
1854  end function sinprod_jac12_square
1855 
1856  function sinprod_jac21_square(eta1, eta2, params)
1857  real(8) :: sinprod_jac21_square
1858  real(8), intent(in) :: eta1
1859  real(8), intent(in) :: eta2
1860  sll_real64, dimension(:), optional, intent(in) :: params
1861  sinprod_jac21_square = 0.5_8*0.1_8*cos(0.5*eta1)*sin(0.5*eta2)
1862  end function sinprod_jac21_square
1863 
1864  function sinprod_jac22_square(eta1, eta2, params)
1865  real(8) :: sinprod_jac22_square
1866  real(8), intent(in) :: eta1
1867  real(8), intent(in) :: eta2
1868  sll_real64, dimension(:), optional, intent(in) :: params
1869  sinprod_jac22_square = 1.0_8 + &
1870  0.5_8*0.1_8*sin(0.5*eta1)*cos(0.5*eta2)
1871  end function sinprod_jac22_square
1872 
1873  ! jacobian ie determinant of jacobian matrix
1874  function sinprod_jac_square(eta1, eta2, params)
1875  real(8) :: sinprod_jac_square
1876  real(8), intent(in) :: eta1
1877  real(8), intent(in) :: eta2
1878  sll_real64, dimension(:), optional, intent(in) :: params
1879  !sll_f_sinprod_jac = 1.0_f64 + 0.2_f64 *sll_p_pi * sin (2*sll_p_pi**(eta1+eta2))
1880  sinprod_jac_square = &
1881  (1.0_8 + 0.5_8*0.1_8*cos(0.5*eta1)*sin(0.5*eta2))* &
1882  (1.0_8 + 0.5_8*0.1_8*sin(0.5*eta1)*cos(0.5*eta2)) - &
1883  0.5_8*0.1_8*sin(0.5*eta1)*cos(0.5*eta2)* &
1884  0.5_8*0.1_8*cos(0.5*eta1)*sin(0.5*eta2)
1885  end function sinprod_jac_square
1886 
1887 #endif
1888 
1889 #if 0
1890  ! What is this???
1891  ! test function
1892  !-------------------
1893  ! direct mapping
1894  function test_x1(eta1, eta2, params)
1895  sll_real64 :: test_x1
1896  sll_real64, intent(in) :: eta1
1897  sll_real64, intent(in) :: eta2
1898  sll_real64, dimension(:), optional, intent(in) :: params
1899  test_x1 = eta1 + 0.1_f64*sin(2.0_f64*sll_p_pi*eta1)
1900  !test_x1 = eta1**2
1901  end function test_x1
1902 
1903  function test_x2(eta1, eta2, params)
1904  sll_real64 :: test_x2
1905  sll_real64, intent(in) :: eta1
1906  sll_real64, intent(in) :: eta2
1907  sll_real64, dimension(:), optional, intent(in) :: params
1908  test_x2 = eta2 + 0.1_f64*sin(2.0_f64*sll_p_pi*eta2)
1909  end function test_x2
1910 
1911  ! inverse mapping
1912  function test_eta1(x1, x2, params)
1913  sll_real64 :: test_eta1
1914  sll_real64, intent(in) :: x1
1915  sll_real64, intent(in) :: x2
1916  sll_real64, dimension(:), optional, intent(in) :: params
1917  test_eta1 = x1/0.1_f64
1918  end function test_eta1
1919 
1920  function test_eta2(x1, x2, params)
1921  sll_real64 :: test_eta2
1922  sll_real64, intent(in) :: x1
1923  sll_real64, intent(in) :: x2
1924  sll_real64, dimension(:), optional, intent(in) :: params
1925  test_eta2 = x2/0.1_f64
1926  end function test_eta2
1927 
1928  ! inverse jacobian matrix
1929  function test_jac11(eta1, eta2, params)
1930  sll_real64 :: test_jac11
1931  sll_real64, intent(in) :: eta1
1932  sll_real64, intent(in) :: eta2
1933  sll_real64, dimension(:), optional, intent(in) :: params
1934  test_jac11 = 1.0_f64/(1.0_f64 + 2.0_f64*sll_p_pi*0.1_f64*cos(2.0_f64*sll_p_pi*eta1))
1935  end function test_jac11
1936 
1937  function test_jac12(eta1, eta2, params)
1938  sll_real64 :: test_jac12
1939  sll_real64, intent(in) :: eta1
1940  sll_real64, intent(in) :: eta2
1941  sll_real64, dimension(:), optional, intent(in) :: params
1942  test_jac12 = 0.0_f64
1943  end function test_jac12
1944 
1945  function test_jac21(eta1, eta2, params)
1946  sll_real64 :: test_jac21
1947  sll_real64, intent(in) :: eta1
1948  sll_real64, intent(in) :: eta2
1949  sll_real64, dimension(:), optional, intent(in) :: params
1950  test_jac21 = 0.0_f64
1951  end function test_jac21
1952 
1953  function test_jac22(eta1, eta2, params)
1954  sll_real64 :: test_jac22
1955  sll_real64, intent(in) :: eta1
1956  sll_real64, intent(in) :: eta2
1957  sll_real64, dimension(:), optional, intent(in) :: params
1958  test_jac22 = 1.0_f64/(1.0_f64 + 2.0_f64*sll_p_pi*0.1_f64*cos(2.0_f64*sll_p_pi*eta2))
1959  end function test_jac22
1960 
1961  ! jacobian ie determinant of jacobian matrix
1962  function test_jac(eta1, eta2, params)
1963  sll_real64 :: test_jac
1964  sll_real64, intent(in) :: eta1
1965  sll_real64, intent(in) :: eta2
1966  sll_real64, dimension(:), optional, intent(in) :: params
1967  test_jac = (1.0_f64 + 2.0_f64*sll_p_pi*0.1_f64*cos(2.0_f64*sll_p_pi*eta1))* &
1968  (1.0_f64 + 2.0_f64*sll_p_pi*0.1_f64*cos(2.0_f64*sll_p_pi*eta2))
1969  !test_jac = 2 * eta1!
1970  end function test_jac
1971 #endif
1972 
1973  ! ***************************************************************************
1974  !
1975  ! Alternative formulation for the polar coordinate transformation:
1976  !
1977  ! X1 = (Rmin + (Rmax-Rmin)*eta1)*cos(2*pi*eta2)
1978  ! X2 = (Rmin + (Rmax-Rmin)*eta1)*sin(2*pi*eta2)
1979  !
1980  ! Where eta1 and eta2 are defined in the interval [0,1]. The 'params' array
1981  ! contains the information (R1, R2). Typically:
1982  !
1983  ! R1 = 0.1
1984  ! R2 = 1.0
1985  !
1986  ! ***************************************************************************
1987 
1989  function sll_f_x1_polar_f(eta1, eta2, params)
1990  sll_real64 :: sll_f_x1_polar_f
1991  sll_real64, intent(in) :: eta1, eta2
1992  sll_real64, dimension(:), intent(in) :: params
1993  sll_real64 :: r1
1994  sll_real64 :: r2
1995 
1996  sll_assert(size(params) >= 2)
1997  r1 = params(1)
1998  r2 = params(2)
1999  sll_f_x1_polar_f = (r1 + (r2 - r1)*eta1)*cos(2.0_f64*sll_p_pi*eta2)
2000  end function sll_f_x1_polar_f
2001 
2003  function sll_f_x2_polar_f(eta1, eta2, params)
2004  sll_real64 :: sll_f_x2_polar_f
2005  sll_real64, intent(in) :: eta1, eta2
2006  sll_real64, dimension(:), intent(in) :: params
2007  sll_real64 :: r1
2008  sll_real64 :: r2
2009 
2010  sll_assert(size(params) >= 2)
2011  r1 = params(1)
2012  r2 = params(2)
2013  sll_f_x2_polar_f = (r1 + (r2 - r1)*eta1)*sin(2.0_f64*sll_p_pi*eta2)
2014  end function sll_f_x2_polar_f
2015 
2017  function sll_f_deriv_x1_polar_f_eta1(eta1, eta2, params)
2018  sll_real64 :: sll_f_deriv_x1_polar_f_eta1
2019  sll_real64, intent(in) :: eta1, eta2
2020  sll_real64, dimension(:), intent(in) :: params
2021  sll_real64 :: r1
2022  sll_real64 :: r2
2023 
2024  sll_assert(size(params) >= 2)
2025  r1 = params(1)
2026  r2 = params(2)
2028  sll_f_deriv_x1_polar_f_eta1 = (r2 - r1)*cos(2.0_f64*sll_p_pi*eta2)
2029  end function sll_f_deriv_x1_polar_f_eta1
2030 
2032  function sll_f_deriv_x1_polar_f_eta2(eta1, eta2, params)
2033  sll_real64 :: sll_f_deriv_x1_polar_f_eta2
2034  sll_real64, intent(in) :: eta1, eta2
2035  sll_real64, dimension(:), intent(in) :: params
2036  sll_real64 :: k
2037  sll_real64 :: r1
2038  sll_real64 :: r2
2039  sll_f_deriv_x1_polar_f_eta2 = eta1 + eta2
2040  sll_assert(size(params) >= 2)
2041  r1 = params(1)
2042  r2 = params(2)
2043  k = 2.0_f64*sll_p_pi
2044  sll_f_deriv_x1_polar_f_eta2 = -(r1 + (r2 - r1)*eta1)*sin(k*eta2)*k
2045  end function sll_f_deriv_x1_polar_f_eta2
2046 
2048  function sll_f_deriv_x2_polar_f_eta1(eta1, eta2, params)
2049  sll_real64 :: sll_f_deriv_x2_polar_f_eta1
2050  sll_real64, intent(in) :: eta1, eta2
2051  sll_real64, dimension(:), intent(in) :: params
2052  sll_real64 :: r1
2053  sll_real64 :: r2
2054 
2055  sll_assert(size(params) >= 2)
2056  r1 = params(1)
2057  r2 = params(2)
2059  sll_f_deriv_x2_polar_f_eta1 = (r2 - r1)*sin(2.0_f64*sll_p_pi*eta2)
2060  end function sll_f_deriv_x2_polar_f_eta1
2061 
2063  function sll_f_deriv_x2_polar_f_eta2(eta1, eta2, params)
2064  sll_real64 :: sll_f_deriv_x2_polar_f_eta2
2065  sll_real64, intent(in) :: eta1, eta2
2066  sll_real64, dimension(:), intent(in) :: params
2067  sll_real64 :: k
2068  sll_real64 :: r1
2069  sll_real64 :: r2
2070  sll_f_deriv_x2_polar_f_eta2 = eta1 + eta2
2071  sll_assert(size(params) >= 2)
2072  r1 = params(1)
2073  r2 = params(2)
2074  k = 2.0_f64*sll_p_pi
2075  sll_f_deriv_x2_polar_f_eta2 = (r1 + (r2 - r1)*eta1)*cos(k*eta2)*k
2076  end function sll_f_deriv_x2_polar_f_eta2
2077 
2079  function sll_f_jacobian_polar_f(eta1, eta2, params) result(jac)
2080  sll_real64 :: jac
2081  sll_real64, intent(in) :: eta1, eta2
2082  sll_real64, dimension(:), intent(in) :: params
2083  sll_real64 :: r1
2084  sll_real64 :: r2
2085 
2086  jac = eta1 + eta2
2087  sll_assert(size(params) >= 2)
2088  r1 = params(1)
2089  r2 = params(2)
2090  jac = 2.0_f64*sll_p_pi*(r1 + (r2 - r1)*eta1)*(r2 - r1)
2091  end function sll_f_jacobian_polar_f
2092 
2093  ! what is the following used for? It is not meant or used for the
2094  ! coordinate transformation class... this one is used in the unit_test_2d.F90
2095  ! file but, what else?
2097  function sll_f_deriv1_jacobian_polar_f(eta1, eta2, params) result(deriv)
2098  sll_real64 :: deriv
2099  sll_real64, intent(in) :: eta1, eta2
2100  sll_real64, dimension(:), intent(in) :: params
2101  sll_real64 :: r1
2102  sll_real64 :: r2
2103 
2104  deriv = eta1 + eta2
2105  sll_assert(size(params) >= 2)
2106  r1 = params(1)
2107  r2 = params(2)
2108  deriv = 2.0_f64*sll_p_pi*(r2 - r1)**2
2109  end function sll_f_deriv1_jacobian_polar_f
2110 
2111 ! why is this here?
2112 !!$ function zero_function(eta1, eta2)
2113 !!$ sll_real64, intent(in) :: eta1, eta2
2114 !!$ sll_real64 :: zero_function
2115 !!$ zero_function = 0.0_f64
2116 !!$ end function zero_function
2117 
2118  !************************************************************************
2119  !
2120  ! 1D maps
2121  !
2122  !************************************************************************
2123 
2124  ! Linear map
2125  !
2126  ! x1 = A + (B-A)*eta1
2127  !
2128  ! The params array is organized as (A,B) with default values:
2129  !
2130  ! A = -1.0
2131  ! B = 1.0
2132 
2134  function linear_map_f(eta, params) result(val)
2135  sll_real64 :: val
2136  sll_real64, intent(in) :: eta
2137  sll_real64, dimension(:), intent(in) :: params
2138  sll_real64 :: a
2139  sll_real64 :: b
2140 
2141  sll_assert(size(params) >= 2)
2142  a = params(1)
2143  b = params(2)
2144  val = (b - a)*eta + a
2145  end function linear_map_f
2146 
2148  function linear_map_jac_f(eta, params) result(val)
2149  sll_real64 :: val
2150  sll_real64, intent(in) :: eta
2151  sll_real64, dimension(:), intent(in) :: params
2152  sll_real64 :: a
2153  sll_real64 :: b
2154 
2155  val = eta
2156  sll_assert(size(params) >= 2)
2157  a = params(1)
2158  b = params(2)
2159  val = (b - a)
2160  end function linear_map_jac_f
2161 
2162 !!$#define A 0.0_f64
2163 !!$#define B 6.2831853071795862_f64
2164 !!$
2165 !!$ function linear_map_poisson_f( eta, params ) result(val)
2166 !!$ sll_real64 :: val
2167 !!$ sll_real64, intent(in) :: eta
2168 !!$ val = (B-A)*eta + A
2169 !!$ end function linear_map_poisson_f
2170 !!$
2171 !!$ function linear_map_poisson_jac_f( eta, params ) result(val)
2172 !!$ sll_real64 :: val
2173 !!$ sll_real64, intent(in) :: eta
2174 !!$ val = (B-A)
2175 !!$ end function linear_map_poisson_jac_f
2176 !!$
2177 !!$#undef A
2178 !!$#undef B
2179 
2180  ! **************************************************************************
2181  !
2182  ! "D-shaped annular geometry";
2183  ! sinusoidal product (see P. Colella et al. JCP 230 (2011) formula
2184  ! (102) p 2968):
2185  ! eta1 = eta1-eta1_min/eta1_max-eta1_min
2186  ! eta2 = eta2-eta2_min/eta2_max-eta2_min
2187  ! x1 = alpha1+(alpha2*(2*eta1-1)+alpha3)*cos(2*pi*eta2+alpha4*sin(2*pi*eta2))
2188  ! x2 = alpha5*(alpha2*(2*eta1-1)+alpha3)*sin(2*pi*eta2)
2189  !
2190  ! By default the values of the alpha parameters are:
2191  ! alpha1 =1.7_f64
2192  ! alpha2 =0.074_f64
2193  ! alpha3 =0.536_f64
2194  ! alpha4 =0.4290421957_f64
2195  ! alpha5 =1.66_f64
2196  !
2197  ! These parameters are stored in the params array as:
2198  ! ( alpha1, alpha2, alpha3, alpha4,alpha5 )
2199  ! Domain: [eta1_min,eta1_max] X [eta2_min,eta2_max]
2200  !
2201  ! **************************************************************************
2202 
2204  function sll_f_sharped_geo_x1(eta1, eta2, params)
2205  sll_real64 :: sll_f_sharped_geo_x1
2206  sll_real64, intent(in) :: eta1
2207  sll_real64, intent(in) :: eta2
2208  sll_real64, dimension(:), intent(in) :: params
2209  sll_real64 :: pi2
2210  sll_real64 :: alpha1
2211  sll_real64 :: alpha2
2212  sll_real64 :: alpha3
2213  sll_real64 :: alpha4
2214  sll_real64 :: alpha5
2215  sll_real64 :: eta1n
2216  sll_real64 :: eta2n
2217  sll_real64 :: eta1_min
2218  sll_real64 :: eta2_min
2219  sll_real64 :: eta1_max
2220  sll_real64 :: eta2_max
2221 
2222  sll_assert(size(params) >= 9)
2223  alpha1 = params(1)
2224  alpha2 = params(2)
2225  alpha3 = params(3)
2226  alpha4 = params(4)
2227  alpha5 = params(5)
2228  eta1_min = params(6)
2229  eta2_min = params(7)
2230  eta1_max = params(8)
2231  eta2_max = params(9)
2232 
2233  pi2 = 2.0_f64*sll_p_pi
2234  eta1n = (eta1 - eta1_min)/(eta1_max - eta1_min)
2235  eta2n = (eta2 - eta2_min)/(eta2_max - eta2_min)
2236  sll_f_sharped_geo_x1 = alpha1 + (alpha2*(2._f64*eta1n - 1._f64) + alpha3)* &
2237  cos(pi2*eta2n + alpha4*sin(pi2*eta2n))
2238  end function sll_f_sharped_geo_x1
2239 
2241  function sll_f_sharped_geo_x2(eta1, eta2, params)
2242  sll_real64 :: sll_f_sharped_geo_x2
2243  sll_real64, intent(in) :: eta1
2244  sll_real64, intent(in) :: eta2
2245  sll_real64, dimension(:), intent(in) :: params
2246  sll_real64 :: pi2
2247  sll_real64 :: alpha1
2248  sll_real64 :: alpha2
2249  sll_real64 :: alpha3
2250  sll_real64 :: alpha4
2251  sll_real64 :: alpha5
2252  sll_real64 :: eta1n
2253  sll_real64 :: eta2n
2254  sll_real64 :: eta1_min
2255  sll_real64 :: eta2_min
2256  sll_real64 :: eta1_max
2257  sll_real64 :: eta2_max
2258 
2259  sll_assert(size(params) >= 9)
2260  alpha1 = params(1)
2261  alpha2 = params(2)
2262  alpha3 = params(3)
2263  alpha4 = params(4)
2264  alpha5 = params(5)
2265  eta1_min = params(6)
2266  eta2_min = params(7)
2267  eta1_max = params(8)
2268  eta2_max = params(9)
2269 
2270  pi2 = 2.0_f64*sll_p_pi
2271  eta1n = (eta1 - eta1_min)/(eta1_max - eta1_min)
2272  eta2n = (eta2 - eta2_min)/(eta2_max - eta2_min)
2273  sll_f_sharped_geo_x2 = alpha5*(alpha2*(2._f64*eta1n - 1._f64) + alpha3)*sin(pi2*eta2n)
2274  end function sll_f_sharped_geo_x2
2275 
2277  function sll_f_sharped_geo_jac11(eta1, eta2, params)
2278  sll_real64 :: sll_f_sharped_geo_jac11
2279  sll_real64, intent(in) :: eta1
2280  sll_real64, intent(in) :: eta2
2281  sll_real64, dimension(:), intent(in) :: params
2282  sll_real64 :: pi2
2283  sll_real64 :: alpha1
2284  sll_real64 :: alpha2
2285  sll_real64 :: alpha3
2286  sll_real64 :: alpha4
2287  sll_real64 :: alpha5
2288  sll_real64 :: eta1n
2289  sll_real64 :: eta2n
2290  sll_real64 :: eta1_min
2291  sll_real64 :: eta2_min
2292  sll_real64 :: eta1_max
2293  sll_real64 :: eta2_max
2294 
2295  sll_assert(size(params) >= 9)
2296  alpha1 = params(1)
2297  alpha2 = params(2)
2298  alpha3 = params(3)
2299  alpha4 = params(4)
2300  alpha5 = params(5)
2301  eta1_min = params(6)
2302  eta2_min = params(7)
2303  eta1_max = params(8)
2304  eta2_max = params(9)
2305 
2306  pi2 = 2.0_f64*sll_p_pi
2307  eta1n = (eta1 - eta1_min)/(eta1_max - eta1_min)
2308  eta2n = (eta2 - eta2_min)/(eta2_max - eta2_min)
2309  sll_f_sharped_geo_jac11 = alpha2*2._f64*cos(pi2*eta2n + alpha4*sin(pi2*eta2n))
2310  sll_f_sharped_geo_jac11 = sll_f_sharped_geo_jac11/(eta1_max - eta1_min)
2311  end function sll_f_sharped_geo_jac11
2312 
2314  function sll_f_sharped_geo_jac12(eta1, eta2, params)
2315  sll_real64 :: sll_f_sharped_geo_jac12
2316  sll_real64, intent(in) :: eta1
2317  sll_real64, intent(in) :: eta2
2318  sll_real64, dimension(:), intent(in) :: params
2319  sll_real64 :: pi2
2320  sll_real64 :: alpha1
2321  sll_real64 :: alpha2
2322  sll_real64 :: alpha3
2323  sll_real64 :: alpha4
2324  sll_real64 :: alpha5
2325  sll_real64 :: eta1n
2326  sll_real64 :: eta2n
2327  sll_real64 :: eta1_min
2328  sll_real64 :: eta2_min
2329  sll_real64 :: eta1_max
2330  sll_real64 :: eta2_max
2331 
2332  sll_assert(size(params) >= 9)
2333  alpha1 = params(1)
2334  alpha2 = params(2)
2335  alpha3 = params(3)
2336  alpha4 = params(4)
2337  alpha5 = params(5)
2338  eta1_min = params(6)
2339  eta2_min = params(7)
2340  eta1_max = params(8)
2341  eta2_max = params(9)
2342 
2343  pi2 = 2.0_f64*sll_p_pi
2344  eta1n = (eta1 - eta1_min)/(eta1_max - eta1_min)
2345  eta2n = (eta2 - eta2_min)/(eta2_max - eta2_min)
2346  sll_f_sharped_geo_jac12 = -(alpha2*(2._f64*eta1n - 1._f64) + alpha3)*(pi2 + pi2*alpha4*cos(pi2*eta2n))* &
2347  sin(pi2*eta2n + alpha4*sin(pi2*eta2n))
2348  sll_f_sharped_geo_jac12 = sll_f_sharped_geo_jac12/(eta2_max - eta2_min)
2349  end function sll_f_sharped_geo_jac12
2350 
2352  function sll_f_sharped_geo_jac21(eta1, eta2, params)
2353  sll_real64 :: sll_f_sharped_geo_jac21
2354  sll_real64, intent(in) :: eta1
2355  sll_real64, intent(in) :: eta2
2356  sll_real64, dimension(:), intent(in) :: params
2357  sll_real64 :: pi2
2358  sll_real64 :: alpha1
2359  sll_real64 :: alpha2
2360  sll_real64 :: alpha3
2361  sll_real64 :: alpha4
2362  sll_real64 :: alpha5
2363  sll_real64 :: eta1n
2364  sll_real64 :: eta2n
2365  sll_real64 :: eta1_min
2366  sll_real64 :: eta2_min
2367  sll_real64 :: eta1_max
2368  sll_real64 :: eta2_max
2369 
2370  sll_assert(size(params) >= 9)
2371  alpha1 = params(1)
2372  alpha2 = params(2)
2373  alpha3 = params(3)
2374  alpha4 = params(4)
2375  alpha5 = params(5)
2376  eta1_min = params(6)
2377  eta2_min = params(7)
2378  eta1_max = params(8)
2379  eta2_max = params(9)
2380  pi2 = 2.0_f64*sll_p_pi
2381  eta1n = (eta1 - eta1_min)/(eta1_max - eta1_min)
2382  eta2n = (eta2 - eta2_min)/(eta2_max - eta2_min)
2383  sll_f_sharped_geo_jac21 = 2._f64*alpha5*alpha2*sin(pi2*eta2n)
2384  sll_f_sharped_geo_jac21 = sll_f_sharped_geo_jac21/(eta1_max - eta1_min)
2385  end function sll_f_sharped_geo_jac21
2386 
2388  function sll_f_sharped_geo_jac22(eta1, eta2, params)
2389  sll_real64 :: sll_f_sharped_geo_jac22
2390  sll_real64, intent(in) :: eta1
2391  sll_real64, intent(in) :: eta2
2392  sll_real64, dimension(:), intent(in) :: params
2393  sll_real64 :: pi2
2394  sll_real64 :: alpha1
2395  sll_real64 :: alpha2
2396  sll_real64 :: alpha3
2397  sll_real64 :: alpha4
2398  sll_real64 :: alpha5
2399  sll_real64 :: eta1n
2400  sll_real64 :: eta2n
2401  sll_real64 :: eta1_min
2402  sll_real64 :: eta2_min
2403  sll_real64 :: eta1_max
2404  sll_real64 :: eta2_max
2405 
2406  sll_assert(size(params) >= 9)
2407  alpha1 = params(1)
2408  alpha2 = params(2)
2409  alpha3 = params(3)
2410  alpha4 = params(4)
2411  alpha5 = params(5)
2412  eta1_min = params(6)
2413  eta2_min = params(7)
2414  eta1_max = params(8)
2415  eta2_max = params(9)
2416  pi2 = 2.0_f64*sll_p_pi
2417  eta1n = (eta1 - eta1_min)/(eta1_max - eta1_min)
2418  eta2n = (eta2 - eta2_min)/(eta2_max - eta2_min)
2419  sll_f_sharped_geo_jac22 = alpha5*(alpha2*(2._f64*eta1n - 1._f64) + alpha3)*cos(pi2*eta2n)*pi2
2420  sll_f_sharped_geo_jac22 = sll_f_sharped_geo_jac22/(eta2_max - eta2_min)
2421  end function sll_f_sharped_geo_jac22
2422 
2424  function d_sharped_geo_jac(eta1, eta2, params)
2425  sll_real64 :: d_sharped_geo_jac
2426  sll_real64, intent(in) :: eta1
2427  sll_real64, intent(in) :: eta2
2428  sll_real64, dimension(:), intent(in) :: params
2429  sll_real64 :: jac_11, jac_12, jac_21, jac_22
2430  sll_real64 :: pi2
2431  sll_real64 :: alpha1
2432  sll_real64 :: alpha2
2433  sll_real64 :: alpha3
2434  sll_real64 :: alpha4
2435  sll_real64 :: alpha5
2436  sll_real64 :: eta1n
2437  sll_real64 :: eta2n
2438  sll_real64 :: eta1_min
2439  sll_real64 :: eta2_min
2440  sll_real64 :: eta1_max
2441  sll_real64 :: eta2_max
2442 
2443  sll_assert(size(params) >= 9)
2444  alpha1 = params(1)
2445  alpha2 = params(2)
2446  alpha3 = params(3)
2447  alpha4 = params(4)
2448  alpha5 = params(5)
2449  eta1_min = params(6)
2450  eta2_min = params(7)
2451  eta1_max = params(8)
2452  eta2_max = params(9)
2453  pi2 = 2.0_f64*sll_p_pi
2454 
2455  eta1n = (eta1 - eta1_min)/(eta1_max - eta1_min)
2456  eta2n = (eta2 - eta2_min)/(eta2_max - eta2_min)
2457 
2458  jac_11 = alpha2*2._f64*cos(pi2*eta2n + alpha4*sin(pi2*eta2n))
2459  jac_12 = -(alpha2*(2._f64*eta1n - 1._f64) + alpha3)*(pi2 + pi2*alpha4*cos(pi2*eta2n))* &
2460  sin(pi2*eta2n + alpha4*sin(pi2*eta2n))
2461  jac_21 = 2._f64*alpha5*alpha2*sin(pi2*eta2n)
2462  jac_22 = alpha5*(alpha2*(2._f64*eta1n - 1._f64) + alpha3)*cos(pi2*eta2n)*pi2
2463  d_sharped_geo_jac = jac_11*jac_22 - jac_12*jac_21
2464 
2465  d_sharped_geo_jac = d_sharped_geo_jac/((eta1_max - eta1_min)*(eta2_max - eta2_min))
2466  end function d_sharped_geo_jac
2467 
2469 
Functions for analytic coordinate transformations.
function, public sll_f_homography_x2(eta1, eta2, params)
direct mapping
function, public sll_f_polygonal_jac21(eta1, eta2, params)
jacobian matrix
function, public sll_f_homography_jac11(eta1, eta2, params)
jacobian maxtrix
function, public sll_f_rubber_sheeting_jac22(eta1, eta2, params)
jacobian maxtrix
function, public sll_f_rubber_sheeting_jac11(eta1, eta2, params)
jacobian maxtrix
function, public sll_f_sinprod_gen_jac22(eta1, eta2, params)
jacobian matrix
function, public sll_f_affine_x1(eta1, eta2, params)
direct mapping
function, public sll_f_polar_shear_jac12(eta1, eta2, params)
jacobian matrix
function, public sll_f_polygonal_jac22(eta1, eta2, params)
jacobian matrix
function, public sll_f_x2_polar_f(eta1, eta2, params)
direct mapping
function, public sll_f_sinprod_jac21(eta1, eta2, params)
jacobian matrix
function, public sll_f_identity_jac21(eta1, eta2, params)
jacobian maxtrix
function polygonal_jac(eta1, eta2, params)
jacobian ie determinant of jacobian matrix not used for the moment (recomputed)
function, public sll_f_sharped_geo_x2(eta1, eta2, params)
direct mapping
function, public sll_f_sinprod_x2(eta1, eta2, params)
direct mapping
function, public sll_f_polygonal_jac11(eta1, eta2, params)
jacobian matrix
function, public sll_f_homography_jac12(eta1, eta2, params)
jacobian maxtrix
function, public sll_f_sinprod_gen_x2(eta1, eta2, params)
direct mapping
function, public sll_f_sinprod_gen_jac21(eta1, eta2, params)
jacobian matrix
function, public sll_f_homography_jac22(eta1, eta2, params)
jacobian maxtrix
function, public sll_f_sinprod_gen_jac11(eta1, eta2, params)
jacobian matrix
function, public sll_f_sinprod_jac(eta1, eta2, params)
jacobian ie determinant of jacobian matrix
function, public sll_f_sinprod_jac22(eta1, eta2, params)
jacobian matrix
function, public sll_f_affine_x2(eta1, eta2, params)
direct mapping
function, public sll_f_rubber_sheeting_x1(eta1, eta2, params)
direct mapping
function, public sll_f_sinprod_x1(eta1, eta2, params)
direct mapping
function, public sll_f_sharped_geo_x1(eta1, eta2, params)
direct mapping
function, public sll_f_polar_jac12(eta1, eta2, params)
jacobian matrix
function, public sll_f_polar_shear_jac11(eta1, eta2, params)
jacobian matrix
function, public sll_f_x1_polar_f(eta1, eta2, params)
direct mapping
function, public sll_f_identity_x1(eta1, eta2, params)
direct mapping
function, public sll_f_sinprod_jac11(eta1, eta2, params)
jacobian matrix
function, public sll_f_affine_jac12(eta1, eta2, params)
jacobian maxtrix
function, public sll_f_identity_jac22(eta1, eta2, params)
jacobian maxtrix
function, public sll_f_polar_jac11(eta1, eta2, params)
jacobian matrix
function, public sll_f_polar_shear_x2(eta1, eta2, params)
direct mapping
function, public sll_f_polar_x1(eta1, eta2, params)
direct mapping
function d_sharped_geo_jac(eta1, eta2, params)
jacobian ie determinant of jacobian matrix
function, public sll_f_sharped_geo_jac12(eta1, eta2, params)
jacobian matrix
function, public sll_f_polar_eta1(x1, x2, params)
inverse mapping
function, public sll_f_identity_jac12(eta1, eta2, params)
jacobian maxtrix
function, public sll_f_rubber_sheeting_jac21(eta1, eta2, params)
jacobian maxtrix
function, public sll_f_homography_x1(eta1, eta2, params)
direct mapping
function, public sll_f_homography_jac21(eta1, eta2, params)
jacobian maxtrix
function, public sll_f_affine_jac21(eta1, eta2, params)
jacobian maxtrix
function, public sll_f_polar_x2(eta1, eta2, params)
direct mapping
function homography_jac(eta1, eta2, params)
jacobian ie determinant of jacobian matrix
function, public sll_f_affine_jac11(eta1, eta2, params)
jacobian maxtrix
function, public sll_f_rubber_sheeting_x2(eta1, eta2, params)
direct mapping
function, public sll_f_polar_shear_x1(eta1, eta2, params)
direct mapping
function, public sll_f_polar_jac21(eta1, eta2, params)
jacobian matrix
function, public sll_f_sinprod_gen_jac12(eta1, eta2, params)
jacobian matrix
function, public sll_f_sinprod_gen_jac(eta1, eta2, params)
jacobian ie determinant of jacobian matrix
function identity_jac(eta1, eta2, params)
jacobian ie determinant of jacobian matrix
function, public sll_f_polar_shear_jac22(eta1, eta2, params)
jacobian matrix
function, public sll_f_sinprod_jac12(eta1, eta2, params)
jacobian matrix
function, public sll_f_affine_jac22(eta1, eta2, params)
jacobian maxtrix
function, public sll_f_sharped_geo_jac22(eta1, eta2, params)
jacobian matrix
function, public sll_f_polygonal_x1(eta1, eta2, params)
direct mapping
function affine_jac(eta1, eta2, params)
jacobian ie determinant of jacobian matrix
function, public sll_f_sharped_geo_jac11(eta1, eta2, params)
jacobian matrix
function rubber_sheeting_jac(eta1, eta2, params)
jacobian ie determinant of jacobian matrix
function, public sll_f_polar_jac22(eta1, eta2, params)
jacobian matrix
function, public sll_f_polar_eta2(x1, x2, params)
inverse mapping
function, public sll_f_sharped_geo_jac21(eta1, eta2, params)
jacobian matrix
function, public sll_f_identity_x2(eta1, eta2, params)
direct mapping
function, public sll_f_identity_jac11(eta1, eta2, params)
jacobian maxtrix
function, public sll_f_sinprod_gen_x1(eta1, eta2, params)
direct mapping
function, public sll_f_polygonal_jac12(eta1, eta2, params)
jacobian matrix
function, public sll_f_polar_shear_jac21(eta1, eta2, params)
jacobian matrix
function polar_jac(eta1, eta2, params)
jacobian ie determinant of jacobian matrix
function, public sll_f_rubber_sheeting_jac12(eta1, eta2, params)
jacobian maxtrix
function, public sll_f_polygonal_x2(eta1, eta2, params)
direct mapping
function polar_shear_jac(eta1, eta2, params)
jacobian ie determinant of jacobian matrix not computed
Fortran module where set some physical and mathematical constants.
real(kind=f64), parameter, public sll_p_pi
    Report Typos and Errors