Changeset 1575 for trunk/LMDZ.COMMON/libf/dynphy_lonlat/grid_atob_m.F90

Timestamp:

Jul 13, 2016, 4:29:03 PM (8 years ago)

Author:

emillour

Message:

All GCMs:
Further adaptations to keep up with changes in LMDZ5 concerning
physics/dynamics separation (up to rev r2575 of LMDZ5)

• dyn3d_common:
• infotrac.F90 : propagate initialisations for INCA (Earth GCM)

• misc:
• wxios.F90: updates to use the XIOS2 library

• dynphy_lonlat:
• grid_atob_m.F90: fix for some zoomed grid interpolation cases

EM

File:

• trunk/LMDZ.COMMON/libf/dynphy_lonlat/grid_atob_m.F90 (modified) (5 diffs)

trunk/LMDZ.COMMON/libf/dynphy_lonlat/grid_atob_m.F90

@@ -6 +6 @@
 
   USE assert_eq_m, ONLY: assert_eq
-  REAL, SAVE :: pi, deg2rad
 
   PRIVATE
@@ -22 +21 @@
 ! Arguments:
   REAL,              INTENT(IN)  :: x_i(:), y_i(:) !-- INPUT  X&Y COOR. (mi)(ni)
-  REAL,              INTENT(IN)  :: x_o(:), y_o(:) !-- OUTPUT X&Y COOR. (mi)(ni)
+  REAL,              INTENT(IN)  :: x_o(:), y_o(:) !-- OUTPUT X&Y COOR. (mo)(no)
   REAL,              INTENT(OUT) :: d_o1(:,:)      !-- OUTPUT FIELD     (mo,no)
   REAL,    OPTIONAL, INTENT(IN)  :: d_i (:,:)      !-- INPUT FIELD      (mi,ni)
-  LOGICAL, OPTIONAL, INTENT(IN)  :: msk (:,:)      !-- MASK             (mi,ni)
+  LOGICAL, OPTIONAL, INTENT(IN)  :: msk (:,:)      !-- MASK             (mo,no)
   REAL,    OPTIONAL, INTENT(OUT) :: d_o2(:,:)      !-- OUTPUT FOR d_i^2 (mo,no)
 !-------------------------------------------------------------------------------
 ! Local variables:
   CHARACTER(LEN=256) :: modname="fine2coarse"
-  INTEGER :: mi, ni, ii, ji, mo, no, io, jo, nr(2), m1, n1, m2, n2, nn
-  INTEGER :: num_tot(SIZE(x_o),SIZE(y_o))
-  LOGICAL :: found(SIZE(x_o),SIZE(y_o)), li
-  LOGICAL :: mask (SIZE(x_i),SIZE(y_i)), lo
-  REAL    :: dist (SIZE(x_o),SIZE(y_o))
-  REAL    :: a(SIZE(x_o)), b(SIZE(x_o)), c(SIZE(y_o)), d(SIZE(y_o)), inc
-  REAL, PARAMETER :: thresh=1.E-5
-!-------------------------------------------------------------------------------
-  mask(:,:)=.TRUE.; IF(PRESENT(msk)) mask(:,:)=msk(:,:)
-  mi=SIZE(x_i); m1=mi; ni=SIZE(y_i); n1=ni
-  mo=SIZE(x_o); m2=mo; no=SIZE(y_o); n2=no
+  INTEGER :: mi, ni, ii, ji, mo, no, io, jo, nr(2), m1,m2, n1,n2, mx,my, nn, i,j
+  LOGICAL :: li, lo, first=.TRUE.
+  REAL    :: inc, cpa, spa, crlo(SIZE(x_i))
+  REAL, SAVE :: pi, hpi
+  INTEGER, DIMENSION(SIZE(x_o),SIZE(y_o)) :: num_tot
+  LOGICAL, DIMENSION(SIZE(x_o),SIZE(y_o)) :: found, mask
+  REAL,    DIMENSION(SIZE(x_i),SIZE(y_i)) :: dist
+  REAL,    DIMENSION(SIZE(x_o))           :: a, b
+  REAL,    DIMENSION(SIZE(y_o))           :: c, d
+  REAL,    PARAMETER :: thresh=1.E-5
+!-------------------------------------------------------------------------------
+  IF(first) THEN; pi=4.0*ATAN(1.0); hpi=pi/2.0; first=.FALSE.; END IF
+  mi=SIZE(x_i); ni=SIZE(y_i); mo=SIZE(x_o); no=SIZE(y_o)
+  m1=m1; m2=mo; mx=mo; IF(PRESENT(msk)) mx=SIZE(msk,1)
+  n1=ni; n2=no; my=no; IF(PRESENT(msk)) my=SIZE(msk,2)
   li=PRESENT(d_i ); IF(li) THEN; m1=SIZE(d_i ,1); n1=SIZE(d_i ,2); END IF
   lo=PRESENT(d_o2); IF(lo) THEN; m2=SIZE(d_o2,1); n2=SIZE(d_o2,2); END IF
-  mi=assert_eq(mi,m1,SIZE(mask,1),TRIM(modname)//" mi")
-  ni=assert_eq(ni,n1,SIZE(mask,2),TRIM(modname)//" ni")
-  mo=assert_eq(mo,m2,SIZE(d_o1,1),TRIM(modname)//" mo")
-  no=assert_eq(no,n2,SIZE(d_o1,2),TRIM(modname)//" no")
+  mi=assert_eq(mi,m1,TRIM(modname)//" mi")
+  ni=assert_eq(ni,n1,TRIM(modname)//" ni")
+  mo=assert_eq(mo,m2,mx,SIZE(d_o1,1),TRIM(modname)//" mo")
+  no=assert_eq(no,n2,my,SIZE(d_o1,2),TRIM(modname)//" no")
+  mask(:,:)=.TRUE.; IF(PRESENT(msk)) mask(:,:)=msk(:,:)
 
 !--- COMPUTE CELLS INTERFACES COORDINATES OF OUTPUT GRID
@@ -67 +71 @@
              (x_i(ii)-a(io)>thresh.OR.x_i(ii)-b(io)<thresh)) CYCLE
           num_tot(io,jo)=num_tot(io,jo)+1
-          IF(mask(ii,ji)) d_o1(io,jo)=d_o1(io,jo)+inc
+          IF(mask(io,jo)) d_o1(io,jo)=d_o1(io,jo)+inc
           IF(.NOT.lo) CYCLE
-          IF(mask(ii,ji)) d_o2(io,jo)=d_o2(io,jo)+inc*inc
+          IF(mask(io,jo)) d_o2(io,jo)=d_o2(io,jo)+inc*inc
         END DO
       END DO
@@ -89 +93 @@
       IF(found(io,jo)) CYCLE
 !      IF(prt_level>=1) PRINT*, "Problem: point out of domain (i,j)=", io,jo
-      CALL dist_sphe(x_o(io),y_o(jo),x_i,y_i,dist(:,:))
+      crlo(:)=COS(x_o(io)-x_i(:))     !--- COS of points 1 and 2 angle
+      cpa=COS(y_o(jo)); spa=SIN(y_o(jo))
+      DO j=1,ni; dist(:,j)=ACOS(spa*SIN(y_i(j))+cpa*COS(y_i(j))*crlo(:)); END DO
       nr=MINLOC(dist(:,:))!; IF(prt_level>=1) PRINT*, "Solution: ", nr
       inc=1.0; IF(li) inc=d_i(nr(1),nr(2))
-      IF(mask(nr(1),nr(2))) d_o1(io,jo)=inc
+      IF(mask(io,jo)) d_o1(io,jo)=inc
     END DO
   END DO
@@ -250 +256 @@
 !-------------------------------------------------------------------------------
 
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE dist_sphe(rf_lon,rf_lat,rlon,rlat,distance)
-!
-!-------------------------------------------------------------------------------
-! Author:  Laurent Li (december 30th 1996).
-! Purpose: Compute min. distance (along big circle) between 2 points in degrees.
-!-------------------------------------------------------------------------------
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-  REAL, INTENT(IN) :: rf_lon, rf_lat  !--- Reference point coordinates (degrees)
-  REAL, INTENT(IN) :: rlon(:), rlat(:)!--- Points longitudes/latitudes (degrees)
-  REAL, INTENT(OUT):: distance(SIZE(rlon),SIZE(rlat)) !--- Distance    (degrees)
-!-------------------------------------------------------------------------------
-! Local variables:
-  LOGICAL, SAVE :: first=.TRUE.
-  REAL    :: pa, pb, cpa, cpab, spa, spab, crlo(SIZE(rlon))
-  INTEGER :: i, j
-!-------------------------------------------------------------------------------
-  IF(first) THEN
-    pi=4.0*ATAN(1.0); deg2rad=pi/180.0; first=.FALSE.
-  END IF
-  crlo(:)=COS((rf_lon-rlon(:))*deg2rad)     !--- COS of points 1 and 2 angle
-  pa=(90.0-rf_lat)*deg2rad                  !--- North Pole - Point 1 distance
-  cpa=COS(pa); spa=SIN(pa)
-  DO j=1,SIZE(rlat)
-    pb=(90.0-rlat(j))*deg2rad               !--- North Pole - Point 2 distance
-    cpab=cpa*COS(pb); spab=spa*SIN(pb)
-    distance(:,j)=ACOS(cpab+spab*crlo(:))/deg2rad
-  END DO
-
-END SUBROUTINE dist_sphe
-!
-!-------------------------------------------------------------------------------
-
 END MODULE grid_atob_m
 !