Vector Rotation - Fortran Code

SUBROUTINE C_FIELD_GEO_TO_MODEL (FIELD_X,
&                                 FIELD_Y)

!-----------------------------------------------------------------------
! purpose: rotate a 2-d vector field from geographic to model coordinates
!          accounting for fact that vectors change orientation
!-----------------------------------------------------------------------
!
! algorithm : project geographic long/lat components onto
!             geographic cartesian x/y/z components, then
!             perform euler angle rotation about three axes obtaining
!             model cartesian x'''/y'''/z''' components, then
!             project model cartesian componets onto
!             model long/lat componets (sorry!)
!
! input :
!
! FIELD_X : x field on model grid but with geographic directions
! FIELD_Y : y field on model grid but with geographic directions
!
! output :
!
! FIELD_X : x field on model grid but with model directions
! FIELD_Y : y field on model grid but with model directions
!
!-----------------------------------------------------------------------

!passed arguments
REAL (TYPE_REAL_8),
&      INTENT (INOUT) ::
&
&      FIELD_X (C_NXL:C_NXH, C_NYL:C_NYH),
&      FIELD_Y (C_NXL:C_NXH, C_NYL:C_NYH)

!local vars
REAL (TYPE_REAL_8) ::
&
&      LON_GEO,
&      LAT_GEO,
&      LON_MOD,
&      LAT_MOD,
&      X_GEO,
&      Y_GEO,
&      Z_GEO,
&      X_MOD,
&      Y_MOD,
&      Z_MOD

!-----------------------------------------------------------------------

!perform at one horizontal grid point at a time
DO J = C_NYL, C_NYH
DO I = C_NXL, C_NXH

!the present vector position in model lat/long coords.
LON_MOD = C_X (I)
LAT_MOD = C_Y (J)

!the present vector position in geogr. lat/long coords.
CALL C_MODEL_TO_GEO (LON_MOD,
&                        LAT_MOD,
&                        LON_GEO,
&                        LAT_GEO)

!convert geogr. long/lat into geogr. cartesian
X_GEO =-FIELD_X (I,J) * C_SIND (LON_GEO) -
&           FIELD_Y (I,J) * C_COSD (LON_GEO) * C_SIND (LAT_GEO)
Y_GEO = FIELD_X (I,J) * C_COSD (LON_GEO) -
&           FIELD_Y (I,J) * C_SIND (LON_GEO) * C_SIND (LAT_GEO)
Z_GEO = FIELD_Y (I,J) * C_COSD (LAT_GEO)

!rotate geogr. cartesian into model cartesian
X_MOD = C_REVERSE (1,1) * X_GEO
&         + C_REVERSE (1,2) * Y_GEO
&         + C_REVERSE (1,3) * Z_GEO
Y_MOD = C_REVERSE (2,1) * X_GEO
&         + C_REVERSE (2,2) * Y_GEO
&         + C_REVERSE (2,3) * Z_GEO
Z_MOD = C_REVERSE (3,1) * X_GEO
&         + C_REVERSE (3,2) * Y_GEO
&         + C_REVERSE (3,3) * Z_GEO

!convert model cartesian into model long/lat
FIELD_X (I,J) =-X_MOD * C_SIND (LON_MOD)
&                 + Y_MOD * C_COSD (LON_MOD)

FIELD_Y (I,J) =-X_MOD * C_COSD (LON_MOD) * C_SIND (LAT_MOD)
&                 - Y_MOD * C_SIND (LON_MOD) * C_SIND (LAT_MOD)
&                 + Z_MOD * C_COSD (LAT_MOD)

ENDDO
ENDDO

!-----------------------------------------------------------------------

END SUBROUTINE C_FIELD_GEO_TO_MODEL