Changeset 2361


Ignore:
Timestamp:
Sep 7, 2015, 5:50:29 PM (9 years ago)
Author:
dcugnet
Message:

In etat0dyn: removed few useless lines: "masque" is always known because etat0dyn is called after etat0phys.
In grid_atob: shape error in routine fine2coarse now fixed: "msk" argument and local variable mask must have the dimensions of the output grid. Working unit of dist_sphe is no longer degree, but radian.

Location:
LMDZ5/trunk/libf/dynlonlat_phylonlat
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/dynlonlat_phylonlat/grid_atob_m.f90

    r2299 r2361  
    66
    77  USE assert_eq_m, ONLY: assert_eq
    8   REAL, SAVE :: pi, deg2rad
    98
    109  PRIVATE
    … …  
    2221! Arguments:
    2322  REAL,              INTENT(IN)  :: x_i(:), y_i(:) !-- INPUT  X&Y COOR. (mi)(ni)
    24   REAL,              INTENT(IN)  :: x_o(:), y_o(:) !-- OUTPUT X&Y COOR. (mi)(ni)
     23  REAL,              INTENT(IN)  :: x_o(:), y_o(:) !-- OUTPUT X&Y COOR. (mo)(no)
    2524  REAL,              INTENT(OUT) :: d_o1(:,:)      !-- OUTPUT FIELD     (mo,no)
    2625  REAL,    OPTIONAL, INTENT(IN)  :: d_i (:,:)      !-- INPUT FIELD      (mi,ni)
    27   LOGICAL, OPTIONAL, INTENT(IN)  :: msk (:,:)      !-- MASK             (mi,ni)
     26  LOGICAL, OPTIONAL, INTENT(IN)  :: msk (:,:)      !-- MASK             (mo,no)
    2827  REAL,    OPTIONAL, INTENT(OUT) :: d_o2(:,:)      !-- OUTPUT FOR d_i^2 (mo,no)
    2928!-------------------------------------------------------------------------------
    3029! Local variables:
    3130  CHARACTER(LEN=256) :: modname="fine2coarse"
    32   INTEGER :: mi, ni, ii, ji, mo, no, io, jo, nr(2), m1, n1, m2, n2, nn
    33   INTEGER :: num_tot(SIZE(x_o),SIZE(y_o))
    34   LOGICAL :: found(SIZE(x_o),SIZE(y_o)), li
    35   LOGICAL :: mask (SIZE(x_i),SIZE(y_i)), lo
    36   REAL    :: dist (SIZE(x_o),SIZE(y_o))
    37   REAL    :: a(SIZE(x_o)), b(SIZE(x_o)), c(SIZE(y_o)), d(SIZE(y_o)), inc
    38   REAL, PARAMETER :: thresh=1.E-5
    39 !-------------------------------------------------------------------------------
    40   mask(:,:)=.TRUE.; IF(PRESENT(msk)) mask(:,:)=msk(:,:)
    41   mi=SIZE(x_i); m1=mi; ni=SIZE(y_i); n1=ni
    42   mo=SIZE(x_o); m2=mo; no=SIZE(y_o); n2=no
     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
     34  INTEGER, DIMENSION(SIZE(x_o),SIZE(y_o)) :: num_tot
     35  LOGICAL, DIMENSION(SIZE(x_o),SIZE(y_o)) :: found, mask
     36  REAL,    DIMENSION(SIZE(x_o),SIZE(y_o)) :: dist
     37  REAL,    DIMENSION(SIZE(x_o))           :: a, b
     38  REAL,    DIMENSION(SIZE(y_o))           :: c, d
     39  REAL,    PARAMETER :: thresh=1.E-5
     40!-------------------------------------------------------------------------------
     41  mi=SIZE(x_i); ni=SIZE(y_i); mo=SIZE(x_o); no=SIZE(y_o)
     42  m1=m1; m2=mo; mx=mo; IF(PRESENT(msk)) mx=SIZE(msk,1)
     43  n1=ni; n2=no; my=no; IF(PRESENT(msk)) my=SIZE(msk,2)
    4344  li=PRESENT(d_i ); IF(li) THEN; m1=SIZE(d_i ,1); n1=SIZE(d_i ,2); END IF
    4445  lo=PRESENT(d_o2); IF(lo) THEN; m2=SIZE(d_o2,1); n2=SIZE(d_o2,2); END IF
    45   mi=assert_eq(mi,m1,SIZE(mask,1),TRIM(modname)//" mi")
    46   ni=assert_eq(ni,n1,SIZE(mask,2),TRIM(modname)//" ni")
    47   mo=assert_eq(mo,m2,SIZE(d_o1,1),TRIM(modname)//" mo")
    48   no=assert_eq(no,n2,SIZE(d_o1,2),TRIM(modname)//" no")
     46  mi=assert_eq(mi,m1,TRIM(modname)//" mi")
     47  ni=assert_eq(ni,n1,TRIM(modname)//" ni")
     48  mo=assert_eq(mo,m2,mx,SIZE(d_o1,1),TRIM(modname)//" mo")
     49  no=assert_eq(no,n2,my,SIZE(d_o1,2),TRIM(modname)//" no")
     50  mask(:,:)=.TRUE.; IF(PRESENT(msk)) mask(:,:)=msk(:,:)
    4951
    5052!--- COMPUTE CELLS INTERFACES COORDINATES OF OUTPUT GRID
    … …  
    6769             (x_i(ii)-a(io)>thresh.OR.x_i(ii)-b(io)<thresh)) CYCLE
    6870          num_tot(io,jo)=num_tot(io,jo)+1
    69           IF(mask(ii,ji)) d_o1(io,jo)=d_o1(io,jo)+inc
     71          IF(mask(io,jo)) d_o1(io,jo)=d_o1(io,jo)+inc
    7072          IF(.NOT.lo) CYCLE
    71           IF(mask(ii,ji)) d_o2(io,jo)=d_o2(io,jo)+inc*inc
     73          IF(mask(io,jo)) d_o2(io,jo)=d_o2(io,jo)+inc*inc
    7274        END DO
    7375      END DO
    … …  
    9294      nr=MINLOC(dist(:,:))!; IF(prt_level>=1) PRINT*, "Solution: ", nr
    9395      inc=1.0; IF(li) inc=d_i(nr(1),nr(2))
    94       IF(mask(nr(1),nr(2))) d_o1(io,jo)=inc
     96      IF(mask(io,jo)) d_o1(io,jo)=inc
    9597    END DO
    9698  END DO
    … …  
    257259!-------------------------------------------------------------------------------
    258260! Author:  Laurent Li (december 30th 1996).
    259 ! Purpose: Compute min. distance (along big circle) between 2 points in degrees.
    260 !-------------------------------------------------------------------------------
    261   IMPLICIT NONE
    262 !-------------------------------------------------------------------------------
    263 ! Arguments:
    264   REAL, INTENT(IN) :: rf_lon, rf_lat  !--- Reference point coordinates (degrees)
    265   REAL, INTENT(IN) :: rlon(:), rlat(:)!--- Points longitudes/latitudes (degrees)
    266   REAL, INTENT(OUT):: distance(SIZE(rlon),SIZE(rlat)) !--- Distance    (degrees)
     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)
    267269!-------------------------------------------------------------------------------
    268270! Local variables:
    269271  LOGICAL, SAVE :: first=.TRUE.
    270   REAL    :: pa, pb, cpa, cpab, spa, spab, crlo(SIZE(rlon))
     272  REAL,    SAVE :: pi, hpi
     273  REAL    :: pa, pb, cpa, spa, crlo(SIZE(rlon))
    271274  INTEGER :: i, j
    272275!-------------------------------------------------------------------------------
    273   IF(first) THEN
    274     pi=4.0*ATAN(1.0); deg2rad=pi/180.0; first=.FALSE.
    275   END IF
    276   crlo(:)=COS((rf_lon-rlon(:))*deg2rad)     !--- COS of points 1 and 2 angle
    277   pa=(90.0-rf_lat)*deg2rad                  !--- North Pole - Point 1 distance
     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
    278279  cpa=COS(pa); spa=SIN(pa)
    279280  DO j=1,SIZE(rlat)
    280     pb=(90.0-rlat(j))*deg2rad               !--- North Pole - Point 2 distance
    281     cpab=cpa*COS(pb); spab=spa*SIN(pb)
    282     distance(:,j)=ACOS(cpab+spab*crlo(:))/deg2rad
     281    pb=hpi-rlat(j)                    !--- North Pole - Point 2 distance
     282    distance(:,j)=ACOS(cpa*COS(pb)+spa*SIN(pb)*crlo(:))
    283283  END DO
    284284

  • LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0dyn_netcdf.F90

    r2336 r2361  
    3636  USE ioipsl,         ONLY: flininfo, flinopen, flinget, flinclo, histclo
    3737  USE assert_eq_m,    ONLY: assert_eq
    38   USE indice_sol_mod, ONLY: epsfra
    3938  IMPLICIT NONE
    4039
    … …  
    8786  INTEGER            :: i, j, l, ji, itau, iday
    8887  REAL               :: xpn, xps, time, phystep
    89   REAL, DIMENSION(iip1,jjp1)       :: psol, masque_tmp
     88  REAL, DIMENSION(iip1,jjp1)       :: psol
    9089  REAL, DIMENSION(iip1,jjp1,llm+1) :: p3d
    9190  REAL, DIMENSION(iip1,jjp1,llm)   :: uvent, t3d, tpot, qsat, qd
    … …  
    101100
    102101  deg2rad = pi/180.0
    103 
    104 ! Compute ground geopotential and possibly the mask.
    105 !*******************************************************************************
    106   masque_tmp(:,:)=masque(:,:)
    107   WRITE(fmt,"(i4,'i1)')")iip1 ; fmt='('//ADJUSTL(fmt)
    108   IF(ALL(masque==-99999.)) THEN                         !--- KEEP NEW MASK
    109     masque=masque_tmp
    110     IF(prt_level>=1) THEN
    111       WRITE(lunout,*)'BUILT MASK :'
    112       WRITE(lunout,fmt) NINT(masque)
    113     END IF
    114     WHERE(   masque(:,:)<EPSFRA) masque(:,:)=0.
    115     WHERE(1.-masque(:,:)<EPSFRA) masque(:,:)=1.
    116   END IF
    117102
    118103! Compute psol AND tsol, knowing phis.
Note: See TracChangeset for help on using the changeset viewer.