Changeset 2533
- Timestamp:
- Jun 1, 2016, 3:51:57 PM (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/branches/testing/libf/dynphy_lonlat/grid_atob_m.f90
r2528 r2533 29 29 ! Local variables: 30 30 CHARACTER(LEN=256) :: modname="fine2coarse" 31 INTEGER :: mi, ni, ii, ji, mo, no, io, jo, nr(2), m1,m2, n1,n2, mx,my, nn, i,j 32 LOGICAL :: li, lo, first=.TRUE. 33 REAL :: inc, cpa, spa, crlo(SIZE(x_i)) 34 REAL, SAVE :: pi, hpi 31 INTEGER :: mi, ni, ii, ji, mo, no, io, jo, nr(2), m1, n1, m2, n2, mx, my, nn 32 LOGICAL :: li, lo 33 REAL :: inc 35 34 INTEGER, DIMENSION(SIZE(x_o),SIZE(y_o)) :: num_tot 36 35 LOGICAL, DIMENSION(SIZE(x_o),SIZE(y_o)) :: found, mask 37 REAL, DIMENSION(SIZE(x_ i),SIZE(y_i)) :: dist36 REAL, DIMENSION(SIZE(x_o),SIZE(y_o)) :: dist 38 37 REAL, DIMENSION(SIZE(x_o)) :: a, b 39 38 REAL, DIMENSION(SIZE(y_o)) :: c, d 40 39 REAL, PARAMETER :: thresh=1.E-5 41 40 !------------------------------------------------------------------------------- 42 IF(first) THEN; pi=4.0*ATAN(1.0); hpi=pi/2.0; first=.FALSE.; END IF43 41 mi=SIZE(x_i); ni=SIZE(y_i); mo=SIZE(x_o); no=SIZE(y_o) 44 42 m1=m1; m2=mo; mx=mo; IF(PRESENT(msk)) mx=SIZE(msk,1) … … 93 91 IF(found(io,jo)) CYCLE 94 92 ! IF(prt_level>=1) PRINT*, "Problem: point out of domain (i,j)=", io,jo 95 crlo(:)=COS(x_o(io)-x_i(:)) !--- COS of points 1 and 2 angle 96 cpa=COS(y_o(jo)); spa=SIN(y_o(jo)) 97 DO j=1,ni; dist(:,j)=ACOS(spa*SIN(y_i(j))+cpa*COS(y_i(j))*crlo(:)); END DO 93 CALL dist_sphe(x_o(io),y_o(jo),x_i,y_i,dist(:,:)) 98 94 nr=MINLOC(dist(:,:))!; IF(prt_level>=1) PRINT*, "Solution: ", nr 99 95 inc=1.0; IF(li) inc=d_i(nr(1),nr(2)) … … 256 252 !------------------------------------------------------------------------------- 257 253 254 255 !------------------------------------------------------------------------------- 256 ! 257 SUBROUTINE dist_sphe(rf_lon,rf_lat,rlon,rlat,distance) 258 ! 259 !------------------------------------------------------------------------------- 260 ! Author: Laurent Li (december 30th 1996). 261 ! Purpose: Compute min. distance (along big circle) between 2 points in radians. 262 !------------------------------------------------------------------------------- 263 IMPLICIT NONE 264 !------------------------------------------------------------------------------- 265 ! Arguments: 266 REAL, INTENT(IN) :: rf_lon, rf_lat !--- Reference point coordinates (radians) 267 REAL, INTENT(IN) :: rlon(:), rlat(:)!--- Points longitudes/latitudes (radians) 268 REAL, INTENT(OUT):: distance(SIZE(rlon),SIZE(rlat)) !--- Distance (radians) 269 !------------------------------------------------------------------------------- 270 ! Local variables: 271 LOGICAL, SAVE :: first=.TRUE. 272 REAL, SAVE :: pi, hpi 273 REAL :: pa, pb, cpa, spa, crlo(SIZE(rlon)) 274 INTEGER :: i, j 275 !------------------------------------------------------------------------------- 276 IF(first) THEN; pi=4.0*ATAN(1.0); hpi=pi/2.0; first=.FALSE.; END IF 277 crlo(:)=COS(rf_lon-rlon(:)) !--- COS of points 1 and 2 angle 278 pa=hpi-rf_lat !--- North Pole - Point 1 distance 279 cpa=COS(pa); spa=SIN(pa) 280 DO j=1,SIZE(rlat) 281 pb=hpi-rlat(j) !--- North Pole - Point 2 distance 282 distance(:,j)=ACOS(cpa*COS(pb)+spa*SIN(pb)*crlo(:)) 283 END DO 284 285 END SUBROUTINE dist_sphe 286 ! 287 !------------------------------------------------------------------------------- 288 258 289 END MODULE grid_atob_m 259 290 !
Note: See TracChangeset
for help on using the changeset viewer.