#!/bin/bash
#===============================================================================
# mapll.sh  -  Conversao entre ponto do modelo Eta e latitude/longitude.
#
#  DIRETO (PE+local -> global -> lat/lon):
#     ./mapll.sh  MYPE  I_local  J_local
#     ./mapll.sh  159   21       25            (o ponto do blow-up)
#     ./mapll.sh  -1    I_global J_global      (modo global direto)
#
#  INVERSO (lat/lon -> MYPE + global (I,J) + local (i,j)):
#     ./mapll.sh  -i  LAT  LON
#     ./mapll.sh  -i  -49.77  -64.36
#
#  Decomposicao override:  INPES=30 JNPES=75 ./mapll.sh ...
#===============================================================================
#--- Grade Eta ams_08km (editar se mudar de configuracao) ---------------------
IM=735 ; JM=1441
TLM0D=-55 ; TPH0D=-19                 # ponto central (lon,lat) geografico
DLMD=.0576923077 ; DPHD=.0538461539   # incrementos da grade rotacionada (graus)
INPES=${INPES:-30} ; JNPES=${JNPES:-60}   # decomposicao (a rodada que caiu = 30x60)
#------------------------------------------------------------------------------
PLOT=0
if [ "$1" = "-p" ]; then PLOT=1 ; shift ; fi   # saida limpa "lat lon" (4 casas) p/ plotar
if [ "$1" = "-i" ] && [ -n "$3" ]; then MODE=INV ; A=$2 ; B=$3
elif [ -n "$3" ];                  then MODE=FWD ; A=$1 ; B=$2 ; C=$3
else
  echo "DIRETO : $0 [-p] MYPE I_local J_local     (ou -1 I_global J_global)"
  echo "INVERSO: $0 [-p] -i LAT LON"
  echo "  -p : saida so com numeros (lat lon | mype gI gJ i j), para pipe/plot"
  exit 1
fi

awk -v IM=$IM -v JM=$JM -v INPES=$INPES -v JNPES=$JNPES -v TLM0D=$TLM0D -v TPH0D=$TPH0D \
    -v DLMD=$DLMD -v DPHD=$DPHD -v MODE=$MODE -v A=$A -v B=$B -v C=$C -v PLOT=$PLOT 'BEGIN{
  d2r = atan2(0,-1)/180.0
  t0  = TPH0D*d2r
  ICHUNK=int(IM/INPES) ; ITAIL=IM%INPES         # 24, 15  (mesma logica do MPPINIT)
  JCHUNK=int(JM/JNPES) ; JTAIL=JM%JNPES         # 24,  1
  WBD = -(IM-1)*DLMD ; SBD = -(JM-1)/2.0*DPHD   # convencao E-grid (corners.f90/const.f)

  if (MODE=="INV") {
    #================= INVERSO: lat/lon -> global -> PE + local ===============
    lat=A ; lon=B
    phi=lat*d2r ; dlon=(lon-TLM0D)*d2r
    Xg=cos(phi)*cos(dlon) ; Yg=cos(phi)*sin(dlon) ; Zg=sin(phi)
    Xr =  cos(t0)*Xg + sin(t0)*Zg               # rotacao inversa (-t0 em torno de Y)
    Yr =  Yg
    Zr = -sin(t0)*Xg + cos(t0)*Zg
    if (Zr> 1) Zr=1 ; if (Zr<-1) Zr=-1
    rlat = atan2(Zr, sqrt(Xr*Xr+Yr*Yr))/d2r
    rlon = atan2(Yr, Xr)/d2r
    gJ = int((rlat-SBD)/DPHD + 0.5) + 1
    gI = int((rlon-WBD-((gJ+1)%2)*DLMD)/(2*DLMD) + 0.5) + 1   # passo 2*DLMD + offset linha par
    if (!PLOT) printf "lat=%.4f lon=%.4f  ->  rotado rlon=%.4f rlat=%.4f\n", lat, lon, rlon, rlat
    if (gI<1 || gI>IM || gJ<1 || gJ>JM) { printf "  GLOBAL: I=%d J=%d  (FORA do dominio %dx%d)\n", gI,gJ,IM,JM ; exit }
    # acha o bloco (PE) dono do ponto e o indice local
    acc=0 ; for(b=1;b<=INPES;b++){ c=ICHUNK ; if(b<=ITAIL)c++ ; if(gI<=acc+c){ipe=b-1 ; li=gI-acc ; break} acc+=c }
    acc=0 ; for(b=1;b<=JNPES;b++){ c=JCHUNK ; if(b<=JTAIL)c++ ; if(gJ<=acc+c){jpe=b-1 ; lj=gJ-acc ; break} acc+=c }
    mype = jpe*INPES + ipe
    if (PLOT) { printf "%d %d %d %d %d\n", mype, gI, gJ, li, lj ; exit }
    printf "  GLOBAL : I=%d  J=%d\n", gI, gJ
    printf "  MYPE   : %d   (ipe=%d, jpe=%d)\n", mype, ipe, jpe
    printf "  LOCAL  : i=%d  j=%d   (no subdominio do PE %d)\n", li, lj, mype
    exit
  }

  #==================== DIRETO: PE+local -> global -> lat/lon =================
  MYPE=A ; ILOC=B ; JLOC=C
  if (MYPE < 0) { gI=ILOC ; gJ=JLOC ; ipe=-1 ; jpe=-1 ; isg=1 ; jsg=1 }
  else {
    if (MYPE >= INPES*JNPES) { print "ERRO: MYPE >= INPES*JNPES=" INPES*JNPES ; exit }
    ipe = MYPE%INPES ; jpe = int(MYPE/INPES)
    isg=1 ; for(b=0;b<ipe;b++){ c=ICHUNK ; if((b+1)<=ITAIL) c++ ; isg+=c }
    jsg=1 ; for(b=0;b<jpe;b++){ c=JCHUNK ; if((b+1)<=JTAIL) c++ ; jsg+=c }
    gI = isg + ILOC - 1 ; gJ = jsg + JLOC - 1
  }
  rlat = SBD + (gJ-1)*DPHD ; rlon = WBD + (gI-1)*2*DLMD + ((gJ+1)%2)*DLMD   # E-grid
  rl=rlon*d2r ; rp=rlat*d2r
  sphi = sin(t0)*cos(rp)*cos(rl) + cos(t0)*sin(rp)
  if (sphi> 1) sphi=1 ; if (sphi<-1) sphi=-1
  glat = atan2(sphi, sqrt(1-sphi*sphi))/d2r
  y = cos(rp)*sin(rl) ; x = cos(t0)*cos(rp)*cos(rl) - sin(t0)*sin(rp)
  glon = (TLM0D*d2r + atan2(y,x))/d2r
  while (glon> 180) glon-=360 ; while (glon<-180) glon+=360
  if (PLOT) { printf "%.4f %.4f\n", glat, glon ; exit }
  if (MYPE>=0) printf "MYPE=%d  (ipe=%d, jpe=%d)   local(i=%d, j=%d)   MY_IS_GLB=%d MY_JS_GLB=%d\n", MYPE,ipe,jpe,ILOC,JLOC,isg,jsg
  printf "  GLOBAL   : I=%d  J=%d      (grade %dx%d)\n", gI, gJ, IM, JM
  printf "  ROTADO   : rlon=%9.4f  rlat=%9.4f\n", rlon, rlat
  printf "  GEOGRAF. : lat=%9.4f  lon=%9.4f\n", glat, glon
}'