hmap_axisNB_eval_all_e Subroutine

PURE SUBROUTINE hmap_axisNB_eval_all_e(xv,q1,q2,dX1_dt,dX2_dt,dX1_dz,dX2_dz, &
                                       Jh,    g_tt,    g_tz,    g_zz,     &
                                       Jh_dq1,g_tt_dq1,g_tz_dq1,g_zz_dq1, &
                                       Jh_dq2,g_tt_dq2,g_tz_dq2,g_zz_dq2, &
                                       g_t1,g_t2,g_z1,g_z2,Gh11,Gh22  )
  ! MODULES
  IMPLICIT NONE
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! INPUT VARIABLES
  TYPE(t_hmap_axisNB_auxvar),INTENT(IN) :: xv    !! precomputed auxiliary variables
  REAL(wp),INTENT(IN)  :: q1,q2       !! solution variables q1,q2
  REAL(wp),INTENT(IN)  :: dX1_dt,dX2_dt  !! theta derivative of solution variables q1,q2
  REAL(wp),INTENT(IN)  :: dX1_dz,dX2_dz  !!  zeta derivative of solution variables q1,q2
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! OUTPUT VARIABLES
  REAL(wp),INTENT(OUT) :: Jh,g_tt,g_tz,g_zz              !! Jac,1/Jac,g_{ab} with a=theta/zeta b=theta/zeta
  REAL(wp),INTENT(OUT) :: Jh_dq1,g_tt_dq1,g_tz_dq1,g_zz_dq1  !! and their variation vs q1
  REAL(wp),INTENT(OUT) :: Jh_dq2,g_tt_dq2,g_tz_dq2,g_zz_dq2  !! and their variation vs q2
  REAL(wp),INTENT(OUT) :: g_t1,g_t2,g_z1,g_z2,Gh11,Gh22  !! dq^{i}/dtheta*G^{i1}, dq^{i}/dtheta*G^{i2}, and dq^{i}/dzeta*G^{i1}, dq^{i}/dzeta*G^{i2} and G^{11},G^{22}
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! LOCAL VARIABLES
  REAL(wp) :: Gh21,Gh31,Gh32,Gh33
  REAL(wp) :: Tq(3)
  !===================================================================================================================================
  ASSOCIATE(  T=>xv%T(:),  N=>xv%N(:),  B=>xv%B(:), Np=>xv%Np(:), Bp=>xv%Bp(:), NxB=>xv%NxB(:),&
             NN=>xv%NN  , BB=>xv%BB  , NB=>xv%NB  ,NpN=>xv%NpN  , &
            NpB=>xv%NpB, BpN=>xv%BpN ,BpB=>xv%BpB   )
  Tq=(T+q1*Np+q2*Bp)
  Gh11=NN
  Gh21=NB
  Gh22=BB
  Gh31 =SUM(N( :)*Tq(:))
  Gh32 =SUM(B( :)*Tq(:))
  Gh33 =SUM(Tq(:)*Tq(:))

  Jh=SUM(Tq*NxB)
  Jh_dq1=SUM(Np*NxB)
  Jh_dq2=SUM(Bp*NxB)

  g_t1 = Gh11 * dX1_dt + Gh21 * dX2_dt
  g_t2 = Gh21 * dX1_dt + Gh22 * dX2_dt
  g_z1 = Gh11 * dX1_dz + Gh21 * dX2_dz + Gh31
  g_z2 = Gh21 * dX1_dz + Gh22 * dX2_dz + Gh32

  g_tt =   dX1_dt *  g_t1         +  dX2_dt *  g_t2
  g_tz =   dX1_dt *  g_z1         +  dX2_dt *  g_z2
  g_zz =   dX1_dz * (g_z1 + Gh31) +  dX2_dz * (g_z2 + Gh32)  + Gh33

  !Gh11/dq1 =0 Gh12/dq1 =0 Gh13/dq1 = Np.N
  !            Gh22/dq1 =0 Gh23/dq1 = Np.B
  !                        Gh33/dq1 = 2*Np.Tq
  !Gh11/dq2 =0 Gh12/dq2 =0 Gh13/dq2 = Bp.N
  !            Gh22/dq2 =0 Gh23/dq2 = Bp.B
  !                        Gh33/dq2 = 2*Bp.Tq
  ! => g_t1 /dq1 =0, g_t1/dq2 =0, g_t2/dq1 =0, g_t2/dq2 =0
  ! => g_z1 /dq1 = Gh31/dq1, g_z1/dq2 =Gh31/dq2, g_z2/dq1 =Gh32/dq1, g_z2/dq2 =Gh32/dq2
  g_tt_dq1 = 0.0_wp
  g_tt_dq2 = 0.0_wp

  g_tz_dq1 = NpN*dX1_dt + NpB*dX2_dt
  g_tz_dq2 = BpN*dX1_dt + BpB*dX2_dt

  g_zz_dq1 = 2.0_wp*(SUM(Np(:)*Tq(:)) + NpN* dX1_dz + NpB* dX2_dz )
  g_zz_dq2 = 2.0_wp*(SUM(Bp(:)*Tq(:)) + BpN* dX1_dz + BpB* dX2_dz )
  END ASSOCIATE
END SUBROUTINE hmap_axisNB_eval_all_e