Changeset 2361 for LMDZ5/trunk/libf/dynlonlat_phylonlat

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:

• grid_atob_m.f90 (modified) (5 diffs)
• phylmd/etat0dyn_netcdf.F90 (modified) (3 diffs)

LMDZ5/trunk/libf/dynlonlat_phylonlat/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, n1, m2, n2, mx, my, nn
+  LOGICAL :: li, lo
+  REAL    :: inc
+  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_o))           :: a, b
+  REAL,    DIMENSION(SIZE(y_o))           :: c, d
+  REAL,    PARAMETER :: thresh=1.E-5
+!-------------------------------------------------------------------------------
+  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 +69 @@
              (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
@@ -92 +94 @@
       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
@@ -257 +259 @@
 !-------------------------------------------------------------------------------
 ! 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)
+! 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    :: pa, pb, cpa, cpab, spa, spab, crlo(SIZE(rlon))
+  REAL,    SAVE :: pi, hpi
+  REAL    :: pa, pb, cpa, spa, 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
+  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=(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
+    pb=hpi-rlat(j)                    !--- North Pole - Point 2 distance
+    distance(:,j)=ACOS(cpa*COS(pb)+spa*SIN(pb)*crlo(:))
   END DO
 

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

@@ -36 +36 @@
   USE ioipsl,         ONLY: flininfo, flinopen, flinget, flinclo, histclo
   USE assert_eq_m,    ONLY: assert_eq
-  USE indice_sol_mod, ONLY: epsfra
   IMPLICIT NONE
 
@@ -87 +86 @@
   INTEGER            :: i, j, l, ji, itau, iday
   REAL               :: xpn, xps, time, phystep
-  REAL, DIMENSION(iip1,jjp1)       :: psol, masque_tmp
+  REAL, DIMENSION(iip1,jjp1)       :: psol
   REAL, DIMENSION(iip1,jjp1,llm+1) :: p3d
   REAL, DIMENSION(iip1,jjp1,llm)   :: uvent, t3d, tpot, qsat, qd
@@ -101 +100 @@
 
   deg2rad = pi/180.0
-
-! Compute ground geopotential and possibly the mask.
-!*******************************************************************************
-  masque_tmp(:,:)=masque(:,:)
-  WRITE(fmt,"(i4,'i1)')")iip1 ; fmt='('//ADJUSTL(fmt)
-  IF(ALL(masque==-99999.)) THEN                         !--- KEEP NEW MASK
-    masque=masque_tmp
-    IF(prt_level>=1) THEN
-      WRITE(lunout,*)'BUILT MASK :'
-      WRITE(lunout,fmt) NINT(masque)
-    END IF
-    WHERE(   masque(:,:)<EPSFRA) masque(:,:)=0.
-    WHERE(1.-masque(:,:)<EPSFRA) masque(:,:)=1.
-  END IF
 
 ! Compute psol AND tsol, knowing phis.