Changeset 2654 in lmdz_wrf

Timestamp:

Jul 2, 2019, 2:08:48 PM (5 years ago)

Author:

lfita

Message:

Adding centered differences version of the vectorial calculus

File:

• trunk/tools/module_scientific.f90 (modified) (7 diffs)

trunk/tools/module_scientific.f90

@@ -19 +19 @@
 ! crossingpoints_polys: Subroutine to provide the crossing points between two polygons
 ! curl2D_1o: Subroutine to compute the first order curl of a 2D vectorial field
+! curl2D_c1o: Subroutine to compute the first order centered curl of a 2D vectorial field
 ! distanceRK: Function to provide the distance between two points
 ! divergence2D_1o: Subroutine to compute the first order divergence of a 2D vectorial field
+! divergence2D_c1o: Subroutine to compute the first order centered divergence of a 2D vectorial field
 ! FindMinimumR_K*: Function returns the location of the minimum in the section between Start and End.
 ! fill3DI_2Dvec: Subroutine to fill a 3D integer matrix from a series of indices from a 2D matrix
 ! fill3DR_2Dvec: Subroutine to fill a 3D float matrix from a series of indices from a 2D matrix
 ! gradient2D_1o: Subroutine to compute the first order gradient of a 2D field
+! gradient2D_c1o: Subroutine to compute the first order centered gradient of a 2D field
 ! grid_within_polygon: Subroutine to determine which grid cells from a matrix lay inside a polygon
 ! grid_spacepercen: Subroutine to compute the space-percentages of a series of grid cells (B) which lay inside another 
@@ -40 +43 @@
 ! intersection_2Dlines: Subroutine to provide the intersection point between two lines on the plane using Cramer's method
 ! join_polygon: Subroutine to join two polygons
+! lap2D_1o: Subroutine to compute the first order laplacian of a 2D vectorial field
+! lap2D_c1o: Subroutine to compute the first order centered laplacian of a 2D vectorial field
 ! look_clockwise_borders: Subroutine to look clock-wise for a next point within a collection of borders 
 !   (limits of a region)
@@ -7607 +7612 @@
   SUBROUTINE gradient2D_1o(dx, dy, var, dsx, dsy, grad)
   ! Subroutine to compute the first order gradient of a 2D field
+  !   FROM: From: https://en.wikipedia.org/wiki/Finite_difference_coefficient
 
     IMPLICIT NONE
@@ -7625 +7631 @@
     fname = 'gradient2D_1o'
 
-    grad = zeroRK
+    grad = fillval64
+
+    DO i=1, dx-1
+      DO j=1, dy-1
+        grad(i,j,:) = (/ (var(i+1,j)-var(i,j))/dsx(i,j), (var(i,j+1)-var(i,j))/dsy(i,j) /)
+      END DO
+    END DO
+
+  END SUBROUTINE gradient2D_1o
+
+  SUBROUTINE gradient2D_c1o(dx, dy, var, dsx, dsy, grad)
+  ! Subroutine to compute the first order centerd gradient of a 2D field
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: dx, dy
+    REAL(r_k), DIMENSION(dx,dy), INTENT(in)              :: var, dsx, dsy
+    REAL(r_k), DIMENSION(dx,dy,2), INTENT(out)           :: grad
+
+! Local 
+    INTEGER                                              :: i, j
+    REAL(r_k)                                            :: d1x, d1y
+
+!!!!!!! Variables
+! dx, dy: shape of the 2D field
+! var: variable
+! dsx, dsy: matrices of distances betweeen grid points along each axis
+! grad: gradient
+
+    fname = 'gradient2D_c1o'
+
+    grad = fillval64
 
     DO i=2, dx-1
       DO j=2, dy-1
-        grad(i,j,:) = (/ (var(i+1,j)-var(i,j))/dsx(i,j), (var(i,j+1)-var(i,j))/dsy(i,j) /)
+        d1x = 0.5*(var(i+1,j)-0.5*var(i-1,j))/dsx(i,j)
+        d1y = 0.5*(var(i,j+1)-0.5*var(i,j-1))/dsy(i,j)
+        grad(i,j,:) = (/ d1x, d1y /)
       END DO
     END DO
 
-  END SUBROUTINE gradient2D_1o
+  END SUBROUTINE gradient2D_c1o
 
   SUBROUTINE divergence2D_1o(dx, dy, dsx, dsy, xvar, yvar, diver)
@@ -7655 +7694 @@
     fname = 'divergence2D_1o'
 
-    diver = zeroRK
+    diver = fillval64
+
+    DO i=1, dx-1
+      DO j=1, dy-1
+        diver(i,j) = (xvar(i+1,j)-xvar(i,j))/dsx(i,j) + (yvar(i,j+1)-yvar(i,j))/dsy(i,j)
+      END DO
+    END DO
+
+  END SUBROUTINE divergence2D_1o
+
+  SUBROUTINE divergence2D_c1o(dx, dy, dsx, dsy, xvar, yvar, diver)
+  ! Subroutine to compute the first order centerd divergence of a 2D vectorial field
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: dx, dy
+    REAL(r_k), DIMENSION(dx,dy), INTENT(in)              :: xvar, yvar, dsx, dsy
+    REAL(r_k), DIMENSION(dx,dy), INTENT(out)             :: diver
+
+! Local 
+    INTEGER                                              :: i, j
+    REAL(r_k)                                            :: d1x, d1y
+
+!!!!!!! Variables
+! dx, dy: shape of the 2D field
+! xvar, yvar: x-y components of vectorial variable
+! dsx, dsy: matrices of distances betweeen grid points along each axis
+! diver: divergence
+
+    fname = 'divergence2D_c1o'
+
+    diver = fillval64
 
     DO i=2, dx-1
       DO j=2, dy-1
-        diver(i,j) = (xvar(i+1,j)-xvar(i,j))/dsx(i,j) + (yvar(i,j+1)-yvar(i,j))/dsy(i,j)
+        d1x = 0.5*(xvar(i+1,j)-0.5*xvar(i-1,j))/dsx(i,j)
+        d1y = 0.5*(yvar(i,j+1)-0.5*yvar(i,j-1))/dsy(i,j)
+        diver(i,j) = d1x + d1y
       END DO
     END DO
 
-  END SUBROUTINE divergence2D_1o
+  END SUBROUTINE divergence2D_c1o
 
   SUBROUTINE curl2D_1o(dx, dy, dsx, dsy, xvar, yvar, curl)
@@ -7685 +7757 @@
     fname = 'curl2D_1o'
 
-    curl = zeroRK
+    curl = fillval64
 
     DO i=1, dx-1
@@ -7695 +7767 @@
   END SUBROUTINE curl2D_1o
 
+  SUBROUTINE curl2D_c1o(dx, dy, dsx, dsy, xvar, yvar, curl)
+  ! Subroutine to compute the first order centered curl of a 2D vectorial field
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: dx, dy
+    REAL(r_k), DIMENSION(dx,dy), INTENT(in)              :: xvar, yvar, dsx, dsy
+    REAL(r_k), DIMENSION(dx,dy), INTENT(out)             :: curl
+
+! Local 
+    INTEGER                                              :: i, j
+    REAL(r_k)                                            :: d1x, d1y
+
+!!!!!!! Variables
+! dx, dy: shape of the 2D field
+! xvar, yvar: x-y components of vectorial variable
+! dsx, dsy: matrices of distances betweeen grid points along each axis
+! curl: curl
+
+    fname = 'curl2D_c1o'
+
+    curl = fillval64
+
+    DO i=1, dx-1
+      DO j=1, dy-1
+        d1x = 0.5*(yvar(i+1,j)-0.5*yvar(i-1,j))/dsx(i,j)
+        d1y = 0.5*(xvar(i,j+1)-0.5*xvar(i,j-1))/dsy(i,j)
+        curl(i,j) = d1x - d1y
+      END DO
+    END DO
+
+  END SUBROUTINE curl2D_c1o
+
+  SUBROUTINE lap2D_1o(dx, dy, dsx, dsy, var, lap)
+  ! Subroutine to compute the first order laplacian of a 2D vectorial field
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: dx, dy
+    REAL(r_k), DIMENSION(dx,dy), INTENT(in)              :: var, dsx, dsy
+    REAL(r_k), DIMENSION(dx,dy), INTENT(out)             :: lap
+
+! Local 
+    INTEGER                                              :: i, j
+    REAL(r_k)                                            :: d2x, d2y
+
+!!!!!!! Variables
+! dx, dy: shape of the 2D field
+! xvar, yvar: x-y components of vectorial variable
+! dsx, dsy: matrices of distances betweeen grid points along each axis
+! lap: laplacian
+
+    fname = 'curl2D_1o'
+
+    lap = fillval64
+
+    DO i=1, dx-2
+      DO j=1, dy-2
+        d2x = (1.*var(i,j)-2.*var(i+1,j)+1.*var(i+2,j))/(dsx(i,j)**2)
+        d2y = (1.*var(i,j)-2.*var(i,j+1)+1.*var(i,j+2))/(dsy(i,j)**2)
+        lap(i,j) = d2x+d2y
+      END DO
+    END DO
+
+  END SUBROUTINE lap2D_1o
+
+
+  SUBROUTINE lap2D_c1o(dx, dy, dsx, dsy, var, lap)
+  ! Subroutine to compute the first order centered laplacian of a 2D vectorial field
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: dx, dy
+    REAL(r_k), DIMENSION(dx,dy), INTENT(in)              :: var, dsx, dsy
+    REAL(r_k), DIMENSION(dx,dy), INTENT(out)             :: lap
+
+! Local 
+    INTEGER                                              :: i, j
+    REAL(r_k)                                            :: d2x, d2y
+
+!!!!!!! Variables
+! dx, dy: shape of the 2D field
+! xvar, yvar: x-y components of vectorial variable
+! dsx, dsy: matrices of distances betweeen grid points along each axis
+! lap: laplacian
+
+    fname = 'curl2D_c1o'
+
+    lap = fillval64
+
+    DO i=2, dx-1
+      DO j=2, dy-1
+        d2x = (1.*var(i-1,j)-2.*var(i,j)+1.*var(i+1,j))/(dsx(i,j)**2)
+        d2y = (1.*var(i,j-1)-2.*var(i,j)+1.*var(i,j+1))/(dsy(i,j)**2)
+        lap(i,j) = d2x+d2y
+      END DO
+    END DO
+
+  END SUBROUTINE lap2D_c1o
+
 END MODULE module_scientific