From 8632e1dff9595f9c67025ee8029d8df1aaadd522 Mon Sep 17 00:00:00 2001
From: Dave Gill <gill@ucar.edu>
Date: Wed, 6 Dec 2017 20:36:14 -0700
Subject: [PATCH] Grid to earth relative is -1 * alpha for Cassini

---
 metgrid/src/rotate_winds_module.F | 17 ++++++-----------
 1 file changed, 6 insertions(+), 11 deletions(-)

diff --git a/metgrid/src/rotate_winds_module.F b/metgrid/src/rotate_winds_module.F
index 848775393..2e9f960af 100644
--- a/metgrid/src/rotate_winds_module.F
+++ b/metgrid/src/rotate_winds_module.F
@@ -161,7 +161,6 @@ subroutine metmap_xform(u, u_mask, v, v_mask, &
             ! Rotate U field
             do j=us2,ue2
                do i=us1,ue1
-
                   diff = idir * (xlon_u(i,j) - proj_stack(current_nest_number)%stdlon)
                   if (diff > 180.) then
                      diff = diff - 360.
@@ -172,8 +171,7 @@ subroutine metmap_xform(u, u_mask, v, v_mask, &
                   ! Calculate the rotation angle, alpha, in radians
                   if (proj_stack(current_nest_number)%code == PROJ_LC) then
                      alpha = diff * proj_stack(current_nest_number)%cone * rad_per_deg * proj_stack(current_nest_number)%hemi 
-                  else if (proj_stack(current_nest_number)%code == PROJ_CASSINI) then
-                     if ( idir == 1 ) then
+                  else if ( (proj_stack(SOURCE_PROJ)%code == PROJ_CASSINI) .and. ( idir == 1 ) ) then
                         call latlon_to_ij ( proj_stack(SOURCE_PROJ) , &
                                             xlat_u(i,j), xlon_u(i,j), x_u, y_u )
                         call ij_to_latlon ( proj_stack(SOURCE_PROJ) , &
@@ -187,10 +185,10 @@ subroutine metmap_xform(u, u_mask, v, v_mask, &
                            diff_lon = diff_lon + 360.
                         end if
                         diff_lat = xlat_u_p1-xlat_u_m1
-                        alpha = atan2(   (-cos(xlat_u(i,j)*rad_per_deg) * diff_lon*rad_per_deg),   &
+                        alpha =-atan2(   (-cos(xlat_u(i,j)*rad_per_deg) * diff_lon*rad_per_deg),   &
                                                         diff_lat*rad_per_deg    &
                                      )
-                     else if ( idir == -1 ) then
+                  else if ( (proj_stack(current_nest_number)%code == PROJ_CASSINI) .and. ( idir == -1 ) ) then
                         if (j == ue2) then
                            diff = xlon_u(i,j)-xlon_u(i,j-1)
                            if (diff > 180.) then
@@ -222,7 +220,6 @@ subroutine metmap_xform(u, u_mask, v, v_mask, &
                                                            (xlat_u(i,j+1)-xlat_u(i,j-1))*rad_per_deg    &   
                                         )   
                         end if
-                     endif
                   else
                      alpha = diff * rad_per_deg * proj_stack(current_nest_number)%hemi 
                   end if
@@ -281,8 +278,7 @@ subroutine metmap_xform(u, u_mask, v, v_mask, &
 
                   if (proj_stack(current_nest_number)%code == PROJ_LC) then
                      alpha = diff * proj_stack(current_nest_number)%cone * rad_per_deg * proj_stack(current_nest_number)%hemi 
-                  else if (proj_stack(current_nest_number)%code == PROJ_CASSINI) then
-                     if ( idir == 1 ) then
+                  else if ( (proj_stack(SOURCE_PROJ)%code == PROJ_CASSINI) .and. ( idir == 1 ) ) then
                         call latlon_to_ij ( proj_stack(SOURCE_PROJ) , &
                                             xlat_v(i,j), xlon_v(i,j), x_v, y_v )
                         call ij_to_latlon ( proj_stack(SOURCE_PROJ) , &
@@ -296,10 +292,10 @@ subroutine metmap_xform(u, u_mask, v, v_mask, &
                            diff_lon = diff_lon + 360.
                         end if
                         diff_lat = xlat_v_p1-xlat_v_m1
-                        alpha = atan2(   (-cos(xlat_v(i,j)*rad_per_deg) * diff_lon*rad_per_deg),   &
+                        alpha =-atan2(   (-cos(xlat_v(i,j)*rad_per_deg) * diff_lon*rad_per_deg),   &
                                                         diff_lat*rad_per_deg    &
                                      )
-                     else if ( idir == -1 ) then
+                  else if ( (proj_stack(current_nest_number)%code == PROJ_CASSINI) .and. ( idir == -1 ) ) then
                         if (j == ve2) then
                            diff = xlon_v(i,j)-xlon_v(i,j-1)
                            if (diff > 180.) then
@@ -331,7 +327,6 @@ subroutine metmap_xform(u, u_mask, v, v_mask, &
                                                            (xlat_v(i,j+1)-xlat_v(i,j-1))*rad_per_deg    &
                                         )
                         end if
-                     end if
                   else
                      alpha = diff * rad_per_deg * proj_stack(current_nest_number)%hemi 
                   end if