Ignore:
Timestamp:
Oct 12, 2016, 2:53:20 PM (8 years ago)
Author:
dcugnet
Message:
  • A (re)startphy.nc file (standard name: "startphy0.nc") can be read by ce0l to get land mask, so mask can be defined (in decreasing priority order) from: 1) "o2a.nc file" if this file is found 2) "startphy0.nc" if this file is found 3) "Relief.nc" otherwise
  • Sub-cell scales parameters for orographic gravity waves can be read from file "oro_params.nc" if the configuration key "read_orop" is TRUE. The effect is to bypass the "grid_noro" routine in ce0l, so that any pre-defined mask (from o2a.nc or startphy0.nc) is then overwritten.
  • The gcm stops if the "limit.nc" records number differs from the current year number of days. A warning is issued in case the gcm calendar does not match the time axis attribute "calendar" (if available) from the "limit.nc" file. This attribute is now added to the "limit.nc" time axis.
  • Few simplifications in grid_noro
  • Few parameters changes in acama_gwd and flott_gwd.
  • Variable d_u can be saved in the outputs.
File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/grid_noro_m.F90

    r2576 r2665  
    66  USE assert_eq_m,       ONLY: assert_eq
    77  PRIVATE
    8   PUBLIC :: grid_noro, grid_noro0
     8  PUBLIC :: grid_noro, grid_noro0, read_noro
    99
    1010
    … …  
    7171! CORRELATIONS OF USN OROGRAPHY GRADIENTS         ! dim (imar+2*iext,jmdp+2)
    7272  REAL, ALLOCATABLE :: zxtzxusn(:,:), zytzyusn(:,:), zxtzyusn(:,:)
    73   REAL, ALLOCATABLE :: mask_tmp(:,:), zmea0(:,:)  ! dim (imar+1,jmar)
    74   REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax)
    75   REAL, ALLOCATABLE :: a(:), b(:)                 ! dim (imax)
    76   REAL, ALLOCATABLE :: c(:), d(:)                 ! dim (jmax)
     73  REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imar+1,jmar)
     74  REAL, ALLOCATABLE :: a(:), b(:)                 ! dim (imar+1)
     75  REAL, ALLOCATABLE :: c(:), d(:)                 ! dim (jmar)
    7776  LOGICAL :: masque_lu
    7877  INTEGER :: i, ii, imdp, imar, iext
    7978  INTEGER :: j, jj, jmdp, jmar, nn
    80   REAL    :: xpi, zdeltax, zlenx, weighx, xincr,  zmeanor0
    81   REAL    :: rad, zdeltay, zleny, weighy, masque, zmeasud0
    82   REAL    :: zbordnor, zmeanor, zstdnor, zsignor, zweinor, zpicnor, zvalnor
    83   REAL    :: zbordsud, zmeasud, zstdsud, zsigsud, zweisud, zpicsud, zvalsud
    84   REAL    :: zbordest, zbordoue, xk, xl, xm, xp, xq, xw
     79  REAL    :: xpi, zdeltax, zlenx, weighx, xincr,  zweinor, xk, xl, xm
     80  REAL    :: rad, zdeltay, zleny, weighy, masque, zweisud, xp, xq, xw
     81
     82
     83
    8584!-------------------------------------------------------------------------------
    8685  imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp")
    … …  
    170169    DO jj = 1, jmar
    171170      DO j = 2,jmdp+1
    172         zlenx  =zleny  *COS(yusn(j))
     171        zlenx=zleny*COS(yusn(j))
    173172        zdeltax=zdeltay*COS(yusn(j))
    174         zbordnor=(xincr+c(jj)-yusn(j))*rad
    175         zbordsud=(xincr-d(jj)+yusn(j))*rad
    176         weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny))
     173        weighy=(xincr+AMIN1(c(jj)-yusn(j),yusn(j)-d(jj)))*rad
     174        weighy=AMAX1(0.,AMIN1(weighy,zleny))
     175
    177176        IF(weighy==0.) CYCLE
    178177        DO i = 2, imdp+2*iext-1
    179           zbordest=(xusn(i)-a(ii)+xincr)*rad*COS(yusn(j))
    180           zbordoue=(b(ii)+xincr-xusn(i))*rad*COS(yusn(j))
    181           weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx))
     178          weighx=(xincr+AMIN1(xusn(i)-a(ii),b(ii)-xusn(i)))*rad*COS(yusn(j))
     179          weighx=AMAX1(0.,AMIN1(weighx,zlenx))
     180
    182181          IF(weighx==0.) CYCLE
    183182          num_tot(ii,jj)=num_tot(ii,jj)+1.0
    … …  
    198197!--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME
    199198  IF(.NOT.masque_lu) THEN
    200     WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
     199    WHERE(weight(:,:)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
    201200  END IF
    202   nn=COUNT(weight(:,1:jmar-1)==0.0)
     201  nn=COUNT(weight(:,:)==0.0)
    203202  IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn
    204203  WHERE(weight(:,:)/=0.0)
    … …  
    224223!--- FIRST FILTER, MOVING AVERAGE OVER 9 POINTS.
    225224!-------------------------------------------------------------------------------
    226   ALLOCATE(zmea0(imar+1,jmar))
    227   zmea0(:,:)=zmea(:,:)                           ! GK211005 (CG) UNSMOOTHED TOPO
     225  zphi(:,:)=zmea(:,:)                           ! GK211005 (CG) UNSMOOTHED TOPO
     226
    228227  CALL MVA9(zmea);  CALL MVA9(zstd);  CALL MVA9(zpic);  CALL MVA9(zval)
    229228  CALL MVA9(zxtzx); CALL MVA9(zxtzy); CALL MVA9(zytzy)
    230229
    231230!--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD. (SURFACE PARAMS MEANINGLESS)
    232   ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0
    233   WHERE(mask>=0.1) mask_tmp = 1.
    234   WHERE(weight(:,:)/=0.0)
    235 !   zphi (:,:)= mask_tmp(:,:)*zmea (:,:) ! GK211005 (CG) not necessarly smoothed
    236     zphi (:,:)= mask_tmp(:,:)*zmea0(:,:)
    237     zmea0(:,:)= mask_tmp(:,:)*zmea0(:,:)
    238     zmea (:,:)= mask_tmp(:,:)*zmea (:,:)
    239     zpic (:,:)= mask_tmp(:,:)*zpic (:,:)
    240     zval (:,:)= mask_tmp(:,:)*zval (:,:)
    241     zstd (:,:)= mask_tmp(:,:)*zstd (:,:)
     231  WHERE(weight(:,:)==0.0.OR.mask<0.1)
     232    zphi(:,:)=0.0; zmea(:,:)=0.0; zpic(:,:)=0.0; zval(:,:)=0.0; zstd(:,:)=0.0
    242233  END WHERE
    243234  DO ii = 1, imar
    244235    DO jj = 1, jmar
    245       IF (weight(ii,jj)/=0.0) THEN
    246       !--- Coefficients K, L et M:
    247         xk=(zxtzx(ii,jj)+zytzy(ii,jj))/2.
    248         xl=(zxtzx(ii,jj)-zytzy(ii,jj))/2.
    249         xm=zxtzy(ii,jj)
    250         xp=xk-SQRT(xl**2+xm**2)
    251         xq=xk+SQRT(xl**2+xm**2)
    252         xw=1.e-8
    253         IF(xp<=xw) xp=0.
    254         IF(xq<=xw) xq=xw
    255         IF(ABS(xm)<=xw) xm=xw*SIGN(1.,xm)
    256       !--- SLOPE
    257         zsig(ii,jj)=SQRT(xq)*mask_tmp(ii,jj)
    258       !---ISOTROPY
    259         zgam(ii,jj)=xp/xq*mask_tmp(ii,jj)
    260       !--- THETA ANGLE
    261         zthe(ii,jj)=57.29577951*ATAN2(xm,xl)/2.*mask_tmp(ii,jj)
    262       END IF
     236      IF(weight(ii,jj)==0.0) CYCLE
     237    !--- Coefficients K, L et M:
     238      xk=(zxtzx(ii,jj)+zytzy(ii,jj))/2.
     239      xl=(zxtzx(ii,jj)-zytzy(ii,jj))/2.
     240      xm=zxtzy(ii,jj)
     241      xp=xk-SQRT(xl**2+xm**2)
     242      xq=xk+SQRT(xl**2+xm**2)
     243      xw=1.e-8
     244      IF(xp<=xw) xp=0.
     245      IF(xq<=xw) xq=xw
     246      IF(ABS(xm)<=xw) xm=xw*SIGN(1.,xm)
     247    !--- SLOPE, ANISOTROPY AND THETA ANGLE
     248      zsig(ii,jj)=SQRT(xq)
     249      zgam(ii,jj)=xp/xq
     250      zthe(ii,jj)=90.*ATAN2(xm,xl)/xpi
    263251    END DO
    264252  END DO
     253  WHERE(weight(:,:)==0.0.OR.mask<0.1)
     254    zsig(:,:)=0.0; zgam(:,:)=0.0; zthe(:,:)=0.0
     255  END WHERE
     256
    265257  WRITE(lunout,*)'  MEAN ORO:' ,MAXVAL(zmea)
    266258  WRITE(lunout,*)'  ST. DEV.:' ,MAXVAL(zstd)
    … …  
    271263  WRITE(lunout,*)'  val:'      ,MAXVAL(zval)
    272264     
    273 !--- Values at poles
    274   zmea0(imar+1,:)=zmea0(1,:)
    275   zmea (imar+1,:)=zmea (1,:)
    276   zphi (imar+1,:)=zphi (1,:)
    277   zpic (imar+1,:)=zpic (1,:)
    278   zval (imar+1,:)=zval (1,:)
    279   zstd (imar+1,:)=zstd (1,:)
    280   zsig (imar+1,:)=zsig (1,:)
    281   zgam (imar+1,:)=zgam (1,:)
    282   zthe (imar+1,:)=zthe (1,:)
    283 
    284   zweinor =SUM(weight(1:imar,   1),DIM=1)
    285   zweisud =SUM(weight(1:imar,jmar),DIM=1)
    286   zmeanor0=SUM(weight(1:imar,   1)*zmea0(1:imar,   1),DIM=1)
    287   zmeasud0=SUM(weight(1:imar,jmar)*zmea0(1:imar,jmar),DIM=1)
    288   zmeanor =SUM(weight(1:imar,   1)*zmea (1:imar,   1),DIM=1)
    289   zmeasud =SUM(weight(1:imar,jmar)*zmea (1:imar,jmar),DIM=1)
    290   zstdnor =SUM(weight(1:imar,   1)*zstd (1:imar,   1),DIM=1)
    291   zstdsud =SUM(weight(1:imar,jmar)*zstd (1:imar,jmar),DIM=1)
    292   zsignor =SUM(weight(1:imar,   1)*zsig (1:imar,   1),DIM=1)
    293   zsigsud =SUM(weight(1:imar,jmar)*zsig (1:imar,jmar),DIM=1)
    294   zpicnor =SUM(weight(1:imar,   1)*zpic (1:imar,   1),DIM=1)
    295   zpicsud =SUM(weight(1:imar,jmar)*zpic (1:imar,jmar),DIM=1)
    296   zvalnor =SUM(weight(1:imar,   1)*zval (1:imar,   1),DIM=1)
    297   zvalsud =SUM(weight(1:imar,jmar)*zval (1:imar,jmar),DIM=1)
    298 
    299   zmea(:,1)=zmeanor /zweinor; zmea(:,jmar)=zmeasud /zweisud
    300 !  zphi(:,1)=zmeanor0/zweinor; zphi(:,jmar)=zmeasud0/zweisud   TO COMMIT
    301   zphi(:,1)=zmeanor /zweinor; zphi(:,jmar)=zmeasud /zweisud
    302   zpic(:,1)=zpicnor /zweinor; zpic(:,jmar)=zpicsud /zweisud
    303   zval(:,1)=zvalnor /zweinor; zval(:,jmar)=zvalsud /zweisud
    304   zstd(:,1)=zstdnor /zweinor; zstd(:,jmar)=zstdsud /zweisud
    305   zsig(:,1)=zsignor /zweinor; zsig(:,jmar)=zsigsud /zweisud
    306   zgam(:,1)=1.;               zgam(:,jmar)=1.
    307   zthe(:,1)=0.;               zthe(:,jmar)=0.
     265!--- Values at redundant longitude
     266  zmea(imar+1,:)=zmea(1,:)
     267  zphi(imar+1,:)=zphi(1,:)
     268  zpic(imar+1,:)=zpic(1,:)
     269  zval(imar+1,:)=zval(1,:)
     270  zstd(imar+1,:)=zstd(1,:)
     271  zsig(imar+1,:)=zsig(1,:)
     272  zgam(imar+1,:)=zgam(1,:)
     273  zthe(imar+1,:)=zthe(1,:)
     274
     275!--- Values at north pole
     276  zweinor  =SUM(weight(1:imar,1))
     277  zmea(:,1)=SUM(weight(1:imar,1)*zmea(1:imar,1))/zweinor
     278  zphi(:,1)=SUM(weight(1:imar,1)*zphi(1:imar,1))/zweinor
     279  zpic(:,1)=SUM(weight(1:imar,1)*zpic(1:imar,1))/zweinor
     280  zval(:,1)=SUM(weight(1:imar,1)*zval(1:imar,1))/zweinor
     281  zstd(:,1)=SUM(weight(1:imar,1)*zstd(1:imar,1))/zweinor
     282  zsig(:,1)=SUM(weight(1:imar,1)*zsig(1:imar,1))/zweinor
     283  zgam(:,1)=1.; zthe(:,1)=0.
     284
     285!--- Values at south pole
     286  zweisud     =SUM(weight(1:imar,jmar),DIM=1)
     287  zmea(:,jmar)=SUM(weight(1:imar,jmar)*zmea(1:imar,jmar))/zweisud
     288  zphi(:,jmar)=SUM(weight(1:imar,jmar)*zphi(1:imar,jmar))/zweisud
     289  zpic(:,jmar)=SUM(weight(1:imar,jmar)*zpic(1:imar,jmar))/zweisud
     290  zval(:,jmar)=SUM(weight(1:imar,jmar)*zval(1:imar,jmar))/zweisud
     291  zstd(:,jmar)=SUM(weight(1:imar,jmar)*zstd(1:imar,jmar))/zweisud
     292  zsig(:,jmar)=SUM(weight(1:imar,jmar)*zsig(1:imar,jmar))/zweisud
     293  zgam(:,jmar)=1.; zthe(:,jmar)=0.
    308294
    309295END SUBROUTINE grid_noro
    … …  
    323309!-------------------------------------------------------------------------------
    324310! Arguments:
    325   REAL, INTENT(IN)   :: xd(:), yd(:) !--- INPUT  COORDINATES     (imdp) (jmdp)
    326   REAL, INTENT(IN)   :: zd(:,:)      !--- INPUT  FIELD           (imdp,jmdp)
    327   REAL, INTENT(IN)   :: x(:), y(:)   !--- OUTPUT COORDINATES     (imar+1) (jmar)
    328   REAL, INTENT(OUT)  :: zphi(:,:)    !--- GEOPOTENTIAL           (imar+1,jmar)
    329   REAL, INTENT(INOUT):: mask(:,:)    !--- MASK                   (imar+1,jmar)
     311  REAL, INTENT(IN)  :: xd(:), yd(:) !--- INPUT  COORDINATES     (imdp) (jmdp)
     312  REAL, INTENT(IN)  :: zd(:,:)      !--- INPUT  FIELD           (imdp,  jmdp)
     313  REAL, INTENT(IN)  :: x(:), y(:)   !--- OUTPUT COORDINATES     (imar+1) (jmar)
     314  REAL, INTENT(OUT) :: zphi(:,:)    !--- GEOPOTENTIAL           (imar+1,jmar)
     315  REAL, INTENT(OUT) :: mask(:,:)    !--- MASK                   (imar+1,jmar)
    330316!-------------------------------------------------------------------------------
    331317! Local variables:
    332318  CHARACTER(LEN=256) :: modname="grid_noro0"
    333319  REAL, ALLOCATABLE :: xusn(:), yusn(:)           ! dim (imdp+2*iext) (jmdp+2)
    334   REAL, ALLOCATABLE :: zusn(:,:)                  ! dim (imdp+2*iext,jmdp+2)
     320  REAL, ALLOCATABLE :: zusn(:,:)                  ! dim (imdp+2*iext,  jmdp+2)
    335321  REAL, ALLOCATABLE :: weight(:,:)                ! dim (imar+1,jmar)
    336   REAL, ALLOCATABLE :: mask_tmp(:,:), zmea(:,:)  ! dim (imar+1,jmar)
    337   REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax)
    338   REAL, ALLOCATABLE :: a(:), b(:)                 ! dim (imax)
    339   REAL, ALLOCATABLE :: c(:), d(:)                 ! dim (jmax)
     322  REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imar+1,jmar)
     323  REAL, ALLOCATABLE :: a(:), b(:)                 ! dim (imar+1)
     324  REAL, ALLOCATABLE :: c(:), d(:)                 ! dim (jmar)
     325
    340326  LOGICAL :: masque_lu
    341327  INTEGER :: i, ii, imdp, imar, iext
    342328  INTEGER :: j, jj, jmdp, jmar, nn
    343   REAL    :: xpi, zlenx, weighx, xincr,  zbordnor, zmeanor, zweinor, zbordest
    344   REAL    :: rad, zleny, weighy, masque, zbordsud, zmeasud, zweisud, zbordoue
     329  REAL    :: xpi, zlenx, zleny, weighx, weighy, xincr, masque, rad
     330
    345331!-------------------------------------------------------------------------------
    346332  imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp")
    … …  
    392378
    393379!--- INITIALIZATIONS:
    394   ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0
    395   ALLOCATE(zmea  (imar+1,jmar)); zmea  (:,:)= 0.0
     380  ALLOCATE(weight(imar+1,jmar)); weight(:,:)=0.0; zphi(:,:)=0.0
    396381
    397382!--- SUMMATION OVER GRIDPOINT AREA
    … …  
    403388    DO jj = 1, jmar
    404389      DO j = 2,jmdp+1
    405         zlenx  =zleny  *COS(yusn(j))
    406         zbordnor=(xincr+c(jj)-yusn(j))*rad
    407         zbordsud=(xincr-d(jj)+yusn(j))*rad
    408         weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny))
    409         IF(weighy/=0) THEN
    410           DO i = 2, imdp+2*iext-1
    411             zbordest=(xusn(i)-a(ii)+xincr)*rad*COS(yusn(j))
    412             zbordoue=(b(ii)+xincr-xusn(i))*rad*COS(yusn(j))
    413             weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx))
    414             IF(weighx/=0)THEN
    415               num_tot(ii,jj)=num_tot(ii,jj)+1.0
    416               IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
    417               weight(ii,jj)=weight(ii,jj)+weighx*weighy
    418               zmea  (ii,jj)=zmea  (ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN
    419             END IF
    420           END DO
    421         END IF
     390        zlenx=zleny*COS(yusn(j))
     391        weighy=(xincr+AMIN1(c(jj)-yusn(j),yusn(j)-d(jj)))*rad
     392        weighy=AMAX1(0.,AMIN1(weighy,zleny))
     393        IF(weighy/=0) CYCLE
     394        DO i = 2, imdp+2*iext-1
     395          weighx=(xincr+AMIN1(xusn(i)-a(ii),b(ii)-xusn(i)))*rad*COS(yusn(j))
     396          weighx=AMAX1(0.,AMIN1(weighx,zlenx))
     397          IF(weighx/=0) CYCLE
     398          num_tot(ii,jj)=num_tot(ii,jj)+1.0
     399          IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
     400          weight(ii,jj)=weight(ii,jj)+weighx*weighy
     401          zphi  (ii,jj)=zphi  (ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN
     402        END DO
    422403      END DO
    423404    END DO
    … …  
    426407!--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME
    427408  IF(.NOT.masque_lu) THEN
    428     WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
     409    WHERE(weight(:,:)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
    429410  END IF
    430   nn=COUNT(weight(:,1:jmar-1)==0.0)
     411  nn=COUNT(weight(:,:)==0.0)
    431412  IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn
    432   WHERE(weight/=0.0) zmea(:,:)=zmea(:,:)/weight(:,:)
     413  WHERE(weight/=0.0) zphi(:,:)=zphi(:,:)/weight(:,:)
    433414
    434415!--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD (<10%: SURF PARAMS MEANINGLESS)
    435   ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0
    436   WHERE(mask>=0.1) mask_tmp = 1.
    437   WHERE(weight(:,:)/=0.0)
    438     zphi(:,:)=mask_tmp(:,:)*zmea(:,:)
    439     zmea(:,:)=mask_tmp(:,:)*zmea(:,:)
    440   END WHERE
     416  WHERE(weight(:,:)==0.0.OR.mask<0.1) zphi(:,:)=0.0
     417  WRITE(lunout,*)'  MEAN ORO:' ,MAXVAL(zphi)
     418
     419!--- Values at redundant longitude and at poles
     420  zphi(imar+1,:)=zphi(1,:)
     421  zphi(:,   1)=SUM(weight(1:imar,   1)*zphi(1:imar,   1))/SUM(weight(1:imar,   1))
     422  zphi(:,jmar)=SUM(weight(1:imar,jmar)*zphi(1:imar,jmar))/SUM(weight(1:imar,jmar))
     423
     424END SUBROUTINE grid_noro0
     425!
     426!-------------------------------------------------------------------------------
     427
     428
     429!-------------------------------------------------------------------------------
     430!
     431SUBROUTINE read_noro(x,y,fname,zphi,zmea,zstd,zsig,zgam,zthe,zpic,zval,mask)
     432!
     433!-------------------------------------------------------------------------------
     434! Purpose: Read parameters usually determined with grid_noro from a file.
     435!===============================================================================
     436  USE netcdf, ONLY: NF90_OPEN,  NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION,        &
     437        NF90_NOERR, NF90_CLOSE, NF90_INQ_VARID, NF90_GET_VAR, NF90_STRERROR,   &
     438        NF90_NOWRITE
     439  IMPLICIT NONE
     440!-------------------------------------------------------------------------------
     441! Arguments:
     442  REAL, INTENT(IN)  :: x(:), y(:)    !--- OUTPUT COORDINATES     (imar+1) (jmar)
     443  CHARACTER(LEN=*), INTENT(IN) :: fname ! PARAMETERS FILE NAME
     444  REAL, INTENT(OUT) :: zphi(:,:)     !--- GEOPOTENTIAL           (imar+1,jmar)
     445  REAL, INTENT(OUT) :: zmea(:,:)     !--- MEAN OROGRAPHY         (imar+1,jmar)
     446  REAL, INTENT(OUT) :: zstd(:,:)     !--- STANDARD DEVIATION     (imar+1,jmar)
     447  REAL, INTENT(OUT) :: zsig(:,:)     !--- SLOPE                  (imar+1,jmar)
     448  REAL, INTENT(OUT) :: zgam(:,:)     !--- ANISOTROPY             (imar+1,jmar)
     449  REAL, INTENT(OUT) :: zthe(:,:)     !--- SMALL AXIS ORIENTATION (imar+1,jmar)
     450  REAL, INTENT(OUT) :: zpic(:,:)     !--- MAXIMUM ALTITUDE       (imar+1,jmar)
     451  REAL, INTENT(OUT) :: zval(:,:)     !--- MINIMUM ALTITUDE       (imar+1,jmar)
     452  REAL, INTENT(OUT) :: mask(:,:)     !--- MASK                   (imar+1,jmar)
     453!-------------------------------------------------------------------------------
     454! Local variables:
     455  CHARACTER(LEN=256) :: modname="read_noro"
     456  INTEGER :: imar, jmar, fid, did, vid
     457  LOGICAL :: masque_lu
     458  REAL :: xpi, d2r
     459!-------------------------------------------------------------------------------
     460  imar=assert_eq([SIZE(x),SIZE(zphi,1),SIZE(zmea,1),SIZE(zstd,1),SIZE(zsig,1), &
     461                          SIZE(zgam,1),SIZE(zthe,1),SIZE(zpic,1),SIZE(zval,1), &
     462                          SIZE(mask,1)],TRIM(modname)//" imar")-1
     463  jmar=assert_eq([SIZE(y),SIZE(zphi,2),SIZE(zmea,2),SIZE(zstd,2),SIZE(zsig,2), &
     464                          SIZE(zgam,2),SIZE(zthe,2),SIZE(zpic,2),SIZE(zval,2), &
     465                          SIZE(mask,2)],TRIM(modname)//" jmar")
     466  xpi=ACOS(-1.0); d2r=xpi/180.
     467  WRITE(lunout,*)"*** Orography parameters at sub-cell scale from file ***"
     468
     469!--- ARE WE USING A READ MASK ?
     470  masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0
     471  WRITE(lunout,*)'Masque lu: ',masque_lu
     472  CALL ncerr(NF90_OPEN(fname,NF90_NOWRITE,fid))
     473  CALL check_dim('x','longitude',x(1:imar))
     474  CALL check_dim('y','latitude' ,y(1:jmar))
     475  IF(.NOT.masque_lu) CALL get_fld('mask',mask)
     476  CALL get_fld('Zphi',zphi)
     477  CALL get_fld('Zmea',zmea)
     478  CALL get_fld('mu'  ,zstd)
     479  CALL get_fld('Zsig',zsig)
     480  CALL get_fld('Zgam',zgam)
     481  CALL get_fld('Zthe',zthe)
     482  zpic=zmea+2*zstd
     483  zval=MAX(0.,zmea-2.*zstd)
     484  CALL ncerr(NF90_CLOSE(fid))
    441485  WRITE(lunout,*)'  MEAN ORO:' ,MAXVAL(zmea)
    442 
    443 !--- Values at poles
    444   zphi(imar+1,:)=zphi(1,:)
    445 
    446   zweinor=SUM(weight(1:imar,   1),DIM=1)
    447   zweisud=SUM(weight(1:imar,jmar),DIM=1)
    448   zmeanor=SUM(weight(1:imar,   1)*zmea(1:imar,   1),DIM=1)
    449   zmeasud=SUM(weight(1:imar,jmar)*zmea(1:imar,jmar),DIM=1)
    450   zphi(:,1)=zmeanor/zweinor; zphi(:,jmar)=zmeasud/zweisud
    451 
    452 END SUBROUTINE grid_noro0
     486  WRITE(lunout,*)'  ST. DEV.:' ,MAXVAL(zstd)
     487  WRITE(lunout,*)'  PENTE:'    ,MAXVAL(zsig)
     488  WRITE(lunout,*)'  ANISOTROP:',MAXVAL(zgam)
     489  WRITE(lunout,*)'  ANGLE:'    ,MINVAL(zthe),MAXVAL(zthe)
     490  WRITE(lunout,*)'  pic:'      ,MAXVAL(zpic)
     491  WRITE(lunout,*)'  val:'      ,MAXVAL(zval)
     492
     493CONTAINS
     494
     495
     496SUBROUTINE get_fld(var,fld)
     497  CHARACTER(LEN=*), INTENT(IN)    :: var
     498  REAL,             INTENT(INOUT) :: fld(:,:)
     499  CALL ncerr(NF90_INQ_VARID(fid,var,vid),var)
     500  CALL ncerr(NF90_GET_VAR(fid,vid,fld(1:imar,:)),var)
     501  fld(imar+1,:)=fld(1,:)
     502END SUBROUTINE get_fld
     503
     504SUBROUTINE check_dim(dimd,nam,dimv)
     505  CHARACTER(LEN=*), INTENT(IN) :: dimd
     506  CHARACTER(LEN=*), INTENT(IN) :: nam
     507  REAL,             INTENT(IN) :: dimv(:)
     508  REAL, ALLOCATABLE :: tmp(:)
     509  INTEGER :: n
     510  CALL ncerr(NF90_INQ_DIMID(fid,dimd,did))
     511  CALL ncerr(NF90_INQUIRE_DIMENSION(fid,did,len=n)); ALLOCATE(tmp(n))
     512  CALL ncerr(NF90_INQ_VARID(fid,dimd,did))
     513  CALL ncerr(NF90_GET_VAR(fid,did,tmp))
     514  IF(MAXVAL(tmp)>xpi) tmp=tmp*d2r
     515  IF(n/=SIZE(dimv).OR.ANY(ABS(tmp-dimv)>1E-6)) THEN
     516    WRITE(lunout,*)'Problem with file "'//TRIM(fname)//'".'
     517    CALL abort_physic(modname,'Grid differs from LMDZ for '//TRIM(nam)//'.',1)
     518  END IF
     519END SUBROUTINE check_dim
     520
     521SUBROUTINE ncerr(ncres,var)
     522  IMPLICIT NONE
     523  INTEGER,          INTENT(IN) :: ncres
     524  CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: var
     525  CHARACTER(LEN=256) :: mess
     526  IF(ncres/=NF90_NOERR) THEN
     527    mess='Problem with file "'//TRIM(fname)//'"'
     528    IF(PRESENT(var)) mess=TRIM(mess)//' and variable "'//TRIM(var)//'"'
     529    WRITE(lunout,*)TRIM(mess)//'.'
     530    CALL abort_physic(modname,NF90_STRERROR(ncres),1)
     531  END IF
     532END SUBROUTINE ncerr
     533
     534END SUBROUTINE read_noro
    453535!
    454536!-------------------------------------------------------------------------------
    … …  
    494576END MODULE grid_noro_m
    495577
     578
Note: See TracChangeset for help on using the changeset viewer.