Changeset 2532 for LMDZ5/trunk/libf/dynphy_lonlat

Timestamp:

Jun 1, 2016, 2:48:21 PM (9 years ago)

Author:

dcugnet

Message:

Bug fix: For high resolutions or particular zoom center locations, an output cell (LMDZ grid) with no matching input cells ("landiceref.nc" file grid) caused a crash.

File:

• LMDZ5/trunk/libf/dynphy_lonlat/grid_atob_m.f90 (modified) (3 diffs)

LMDZ5/trunk/libf/dynphy_lonlat/grid_atob_m.f90

@@ -29 +29 @@
 ! Local variables:
   CHARACTER(LEN=256) :: modname="fine2coarse"
-  INTEGER :: mi, ni, ii, ji, mo, no, io, jo, nr(2), m1, n1, m2, n2, mx, my, nn
-  LOGICAL :: li, lo
-  REAL    :: inc
+  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_o),SIZE(y_o)) :: dist
+  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)
@@ -91 +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))
@@ -252 +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 radians.
-!-------------------------------------------------------------------------------
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-  REAL, INTENT(IN) :: rf_lon, rf_lat  !--- Reference point coordinates (radians)
-  REAL, INTENT(IN) :: rlon(:), rlat(:)!--- Points longitudes/latitudes (radians)
-  REAL, INTENT(OUT):: distance(SIZE(rlon),SIZE(rlat)) !--- Distance    (radians)
-!-------------------------------------------------------------------------------
-! Local variables:
-  LOGICAL, SAVE :: first=.TRUE.
-  REAL,    SAVE :: pi, hpi
-  REAL    :: pa, pb, cpa, spa, crlo(SIZE(rlon))
-  INTEGER :: i, j
-!-------------------------------------------------------------------------------
-  IF(first) THEN; pi=4.0*ATAN(1.0); hpi=pi/2.0; first=.FALSE.; END IF
-  crlo(:)=COS(rf_lon-rlon(:))         !--- COS of points 1 and 2 angle
-  pa=hpi-rf_lat                       !--- North Pole - Point 1 distance
-  cpa=COS(pa); spa=SIN(pa)
-  DO j=1,SIZE(rlat)
-    pb=hpi-rlat(j)                    !--- North Pole - Point 2 distance
-    distance(:,j)=ACOS(cpa*COS(pb)+spa*SIN(pb)*crlo(:))
-  END DO
-
-END SUBROUTINE dist_sphe
-!
-!-------------------------------------------------------------------------------
-
 END MODULE grid_atob_m
 !