move compute kernels from physics.f90 to where they belong

src/diagnostics/wind_from_lonlat.F90

@@ -4,11 +4,22 @@ MODULE wind_from_lonlat_mod
   IMPLICIT NONE
   PRIVATE
 
-  PUBLIC :: ulonlat2un, compute_wind2d_perp_from_lonlat_compound, &
+  PUBLIC :: add_lonlat_wind_onto_edge, &
+       ulonlat2un, compute_wind2d_perp_from_lonlat_compound, &
        compute_wind_centered_from_lonlat_compound, compute_wind_perp_from_lonlat_compound       
        
 CONTAINS
 
+  SUBROUTINE add_lonlat_wind_onto_edge(factor, dulon, dulat, ue)
+    USE wind_on_edge_mod, ONLY : add_centered_wind_onto_edge
+    REAL(rstd), INTENT(IN)    :: factor
+    REAL(rstd), INTENT(IN)    :: dulon(:,:), dulat(:,:)
+    REAL(rstd), INTENT(INOUT) :: ue(:,:)
+    REAL(rstd) :: duc(SIZE(dulon,1), SIZE(dulon, 2), 3)
+    CALL compute_wind_centered_from_lonlat_compound(dulon, dulat, duc)
+    CALL add_centered_wind_onto_edge(factor, duc, ue)
+  END SUBROUTINE add_lonlat_wind_onto_edge
+
   SUBROUTINE ulonlat2un(f_ulon, f_ulat,f_u)
     TYPE(t_field), POINTER :: f_ulon(:), f_ulat(:) ! IN : velocity reconstructed at hexagons
     TYPE(t_field), POINTER :: f_u(:)               ! OUT  : normal velocity components on edges

src/diagnostics/wind_on_edge.F90

@@ -4,10 +4,31 @@ MODULE wind_on_edge_mod
   IMPLICIT NONE
   PRIVATE
   
-  PUBLIC :: compute_wind_on_edge
+  PUBLIC :: compute_wind_on_edge, add_centered_wind_onto_edge
   
 CONTAINS
   
+  SUBROUTINE add_centered_wind_onto_edge(factor, uc, ue)
+    REAL(rstd), INTENT(IN)    :: factor
+    REAL(rstd), INTENT(IN)    :: uc(iim*jjm,llm,3)
+    REAL(rstd), INTENT(INOUT) :: ue(3*iim*jjm,llm)
+    INTEGER :: ij,l
+    REAL(rstd) :: due, fac
+    fac = 0.5*factor
+    DO l=ll_begin,ll_end
+       DO ij=ij_begin, ij_end
+          due = sum( (uc(ij,l,:) + uc(ij+t_right,l,:))*ep_e(ij+u_right,:) )
+          ue(ij+u_right,l) = ue(ij+u_right,l) + fac*due
+          
+          due = sum( (uc(ij,l,:) + uc(ij+t_lup,l,:))*ep_e(ij+u_lup,:) )
+          ue(ij+u_lup,l)=ue(ij+u_lup,l) + fac*due
+          
+          due = sum( (uc(ij,l,:) + uc(ij+t_ldown,l,:))*ep_e(ij+u_ldown,:) )
+          ue(ij+u_ldown,l)=ue(ij+u_ldown,l) + fac*due
+       ENDDO
+    ENDDO
+  END SUBROUTINE add_centered_wind_onto_edge
+
   SUBROUTINE compute_wind_on_edge(ue,uedge)
     USE grid_param
     REAL(rstd) :: ue(3*iim*jjm,llm)
@@ -87,5 +108,6 @@ CONTAINS
     END DO
     
   END SUBROUTINE compute_tangential_compound
-  
+
+
 END MODULE wind_on_edge_mod

src/physics/physics.f90

@@ -219,6 +219,7 @@ CONTAINS
     USE physics_dcmip2016_mod, ONLY : full_physics_dcmip2016 => full_physics
     USE etat0_venus_mod, ONLY : full_physics_venus=>full_physics
     USE theta2theta_rhodz_mod
+    USE wind_from_lonlat_mod, ONLY : add_lonlat_wind_onto_edge
     USE mpipara
     USE checksum_mod
     TYPE(t_field),POINTER :: f_phis(:)
@@ -302,7 +303,7 @@ CONTAINS
        ue=f_ue(ind)
        dulon=f_dulon(ind)
        dulat=f_dulat(ind)
-       CALL compute_update_velocity(dulon, dulat, ue)
+       CALL add_lonlat_wind_onto_edge(physics_inout%dt_phys, dulon, dulat, ue)
     END DO
 
   END SUBROUTINE physics_column
@@ -360,32 +361,4 @@ CONTAINS
 
   END SUBROUTINE unpack_physics
 
-  SUBROUTINE compute_update_velocity(dulon, dulat, ue)
-    USE wind_from_lonlat_mod
-    REAL(rstd) :: dulon(iim*jjm,llm)
-    REAL(rstd) :: dulat(iim*jjm,llm)
-    REAL(rstd) :: ue(3*iim*jjm,llm)
-    REAL(rstd) :: duc(iim*jjm,llm,3)
-    REAL(rstd) :: dt2, due
-    INTEGER :: i,j,ij,l
-    ! Reconstruct wind tendencies at edges and add to normal wind
-    CALL compute_wind_centered_from_lonlat_compound(dulon,dulat,duc)
-    dt2=.5*physics_inout%dt_phys
-    DO l=ll_begin,ll_end
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          ij=(j-1)*iim+i
-          due = sum( (duc(ij,l,:) + duc(ij+t_right,l,:))*ep_e(ij+u_right,:) )
-          ue(ij+u_right,l) = ue(ij+u_right,l) + dt2*due
-
-          due = sum( (duc(ij,l,:) + duc(ij+t_lup,l,:))*ep_e(ij+u_lup,:) )
-          ue(ij+u_lup,l)=ue(ij+u_lup,l) + dt2*due
-
-          due = sum( (duc(ij,l,:) + duc(ij+t_ldown,l,:))*ep_e(ij+u_ldown,:) )
-          ue(ij+u_ldown,l)=ue(ij+u_ldown,l) + dt2*due
-        ENDDO
-      ENDDO
-    ENDDO
-  END SUBROUTINE compute_update_velocity
-
 END MODULE physics_mod