Changeset 2576


Ignore:
Timestamp:
Jun 20, 2016, 7:29:33 PM (8 years ago)
Author:
dcugnet
Message:

Use "traditional" sst/sic files without additional records from contiguous years,
but check for previous and next files. Example for Amip Sea Surface Temperature:

  • amipbc_sst_1x1.nc must be present.
  • if amipbc_sst_1x1_m.nc is present (previous year), its last record is used

for interpolation ; otherwise, the last record of amipbc_sst_1x1.nc is used.

  • if amipbc_sst_1x1_p.nc is present (next year), its first record is used

for interpolation ; otherwise, the first record of amipbc_sst_1x1.nc is used.
If both contiguous fiels are missing, the interpolation then implies periodicity.
Note: grid_noro0 has been moved from limit_netcd to grid_noro_m where it belongs.

Location:
LMDZ5/trunk/libf
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/dynphy_lonlat/phylmd/limit_netcdf.F90

    r2540 r2576  
    117117    ALLOCATE(pctsrf(klon,nbsrf))
    118118    CALL start_init_subsurf(.FALSE.)
     119  !--- TO MATCH EXACTLY WHAT WOULD BE DONE IN etat0phys_netcdf
     120    WHERE(   masque(:,:)<EPSFRA) masque(:,:)=0.
     121    WHERE(1.-masque(:,:)<EPSFRA) masque(:,:)=1.
    119122  END IF
    120123
    … …  
    336339  REAL, ALLOCATABLE :: champan(:,:,:)
    337340!--- input files
     341  CHARACTER(LEN=20) :: fnam_m, fnam_p     ! previous/next files names
    338342  CHARACTER(LEN=20) :: cal_in             ! calendar
    339343  CHARACTER(LEN=20) :: unit_sic           ! attribute unit in sea-ice file
    … …  
    345349  LOGICAL           :: extrp              ! flag for extrapolation
    346350  REAL              :: chmin, chmax
    347   INTEGER ierr
     351  INTEGER ierr, idx
    348352  integer n_extrap ! number of extrapolated points
    349353  logical skip
    … …  
    360364  END SELECT
    361365  extrp=.FALSE.; IF(PRESENT(flag).AND.mode=='SST') extrp=flag
     366  idx=INDEX(fnam,'.nc')-1
    362367
    363368!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
    … …  
    393398!--- Time (variable is not needed - it is rebuilt - but calendar is)
    394399  CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(3), name=dnam, len=lmdep), fnam)
    395   ALLOCATE(timeyear(MAX(lmdep,14)))
     400  ALLOCATE(timeyear(lmdep+2))
    396401  CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
    397402  cal_in=' '
    … …  
    412417  IF(lmdep==12) THEN
    413418    timeyear=mid_month(anneeref, cal_in, ndays_in)
    414     CALL msg(0,'Monthly periodic input file (perpetual run).')
    415   ELSE IF(lmdep==14) THEN
    416     timeyear=mid_month(anneeref, cal_in, ndays_in)
    417     CALL msg(0,'Monthly 14-records input file (interannual run).')
     419    CALL msg(0,'Monthly input file(s) for '//TRIM(title)//'.')
    418420  ELSE IF(lmdep==ndays_in) THEN
    419     timeyear=[(REAL(k)-0.5,k=1,ndays_in)]
     421    timeyear=[(REAL(k)-0.5,k=0,ndays_in+1)]
    420422    CALL msg(0,'Daily input file (no time interpolation).')
    421423  ELSE
    422424    WRITE(mess,'(a,i3,a,i3,a)')'Mismatching input file: found',lmdep,        &
    423       ' records, 12/14/',ndays_in,' (periodic/interannual/daily) needed'
     425      ' records, 12/',ndays_in,' (monthly/daily needed).'
    424426    CALL abort_physic('mid_months',TRIM(mess),1)
    425427  END IF
    426428
    427429!--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
    428   ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, MAX(lmdep,14)))
     430  ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep+2))
    429431  IF(extrp) ALLOCATE(work(imdep, jmdep))
    430432  CALL msg(5,'')
    … …  
    446448      WHERE(NINT(mask)/=1) champint=0.001
    447449    END IF
    448     !--- Special case for periodic input file: index shifted
    449     ll = l; IF(lmdep==12) ll = l + 1
    450     champtime(:, :, ll)=champint
     450    champtime(:, :, l+1)=champint
    451451  END DO
    452   IF(lmdep==12) THEN
    453     champtime(:,:, 1)=champtime(:,:,13)
    454     champtime(:,:,14)=champtime(:,:, 2)
    455   END IF
    456452  CALL ncerr(NF90_CLOSE(ncid), fnam)
    457453
     454!--- FIRST RECORD: LAST ONE OF PREVIOUS YEAR (CURRENT YEAR IF UNAVAILABLE)
     455  fnam_m=fnam(1:idx)//'_m.nc'
     456  IF(NF90_OPEN(fnam_m,NF90_NOWRITE,ncid)==NF90_NOERR) THEN
     457    CALL msg(0,'Reading previous year file ("'//TRIM(fnam_m)//'") last record for '//TRIM(title))
     458    CALL ncerr(NF90_INQ_VARID(ncid, varname, varid),fnam_m)
     459    CALL ncerr(NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids),fnam_m)
     460    CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(3), len=l), fnam_m)
     461    CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,l],[imdep,jmdep,1]),fnam_m)
     462    CALL ncerr(NF90_CLOSE(ncid), fnam_m)
     463    CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.)
     464    IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
     465    IF(mode=='RUG') champ=LOG(champ)
     466    CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint)
     467    IF(mode=='RUG') THEN
     468      champint=EXP(champint)
     469      WHERE(NINT(mask)/=1) champint=0.001
     470    END IF
     471    champtime(:, :, 1)=champint
     472  ELSE
     473    CALL msg(0,'Using current year file ("'//TRIM(fnam)//'") last record for '//TRIM(title))
     474    champtime(:, :, 1)=champtime(:, :, lmdep+1)
     475  END IF
     476
     477!--- LAST RECORD: FIRST ONE OF NEXT YEAR (CURRENT YEAR IF UNAVAILABLE)
     478  fnam_p=fnam(1:idx)//'_p.nc'
     479  IF(NF90_OPEN(fnam_p,NF90_NOWRITE,ncid)==NF90_NOERR) THEN
     480    CALL msg(0,'Reading previous year file ("'//TRIM(fnam_p)//'") first record for '//TRIM(title))
     481    CALL ncerr(NF90_INQ_VARID(ncid, varname, varid),fnam_p)
     482    CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,1],[imdep,jmdep,1]),fnam_p)
     483    CALL ncerr(NF90_CLOSE(ncid), fnam_p)
     484    CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.)
     485    IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
     486    IF(mode=='RUG') champ=LOG(champ)
     487    CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint)
     488    IF(mode=='RUG') THEN
     489      champint=EXP(champint)
     490      WHERE(NINT(mask)/=1) champint=0.001
     491    END IF
     492    champtime(:, :, lmdep+2)=champint
     493  ELSE
     494    CALL msg(0,'Using current year file ("'//TRIM(fnam)//'") first record for '//TRIM(title))
     495    champtime(:, :, lmdep+2)=champtime(:, :, 2)
     496  END IF
    458497  DEALLOCATE(dlon_ini, dlat_ini, dlon, dlat, champ)
    459498  IF(extrp) DEALLOCATE(work)
    … …  
    477516     END IF
    478517  END IF
    479   ALLOCATE(yder(MAX(lmdep,14)), champan(iip1, jjp1, ndays))
     518  ALLOCATE(yder(lmdep+2), champan(iip1, jjp1, ndays))
    480519  IF(lmdep==ndays_in) THEN
    481520     champan(1:iim,:,:)=champtime
    … …  
    546585!
    547586!-------------------------------------------------------------------------------
     587  USE grid_noro_m, ONLY: grid_noro0
    548588  IMPLICIT NONE
    549589!===============================================================================
    … …  
    599639
    600640END SUBROUTINE start_init_orog0
    601 !
    602 !-------------------------------------------------------------------------------
    603 
    604 
    605 !-------------------------------------------------------------------------------
    606 !
    607 SUBROUTINE grid_noro0(xd,yd,zd,x,y,zphi,mask)
    608 !
    609 !===============================================================================
    610 ! Purpose: Extracted from grid_noro to provide geopotential height for dynamics
    611 !          without any call to physics subroutines.
    612 !===============================================================================
    613   IMPLICIT NONE
    614 !-------------------------------------------------------------------------------
    615 ! Arguments:
    616   REAL, INTENT(IN)   :: xd(:), yd(:) !--- INPUT  COORDINATES     (imdp) (jmdp)
    617   REAL, INTENT(IN)   :: zd(:,:)      !--- INPUT  FIELD           (imdp,jmdp)
    618   REAL, INTENT(IN)   :: x(:), y(:)   !--- OUTPUT COORDINATES     (imar+1) (jmar)
    619   REAL, INTENT(OUT)  :: zphi(:,:)    !--- GEOPOTENTIAL           (imar+1,jmar)
    620   REAL, INTENT(INOUT):: mask(:,:)    !--- MASK                   (imar+1,jmar)
    621 !-------------------------------------------------------------------------------
    622 ! Local variables:
    623   CHARACTER(LEN=256) :: modname="grid_noro0"
    624   REAL, ALLOCATABLE :: xusn(:), yusn(:)           ! dim (imdp+2*iext) (jmdp+2)
    625   REAL, ALLOCATABLE :: zusn(:,:)                  ! dim (imdp+2*iext,jmdp+2)
    626   REAL, ALLOCATABLE :: weight(:,:)                ! dim (imar+1,jmar)
    627   REAL, ALLOCATABLE :: mask_tmp(:,:), zmea(:,:)   ! dim (imar+1,jmar)
    628   REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax)
    629   REAL, ALLOCATABLE :: a(:), b(:)                 ! dim (imax)
    630   REAL, ALLOCATABLE :: c(:), d(:)                 ! dim (jmax)
    631   LOGICAL :: masque_lu
    632   INTEGER :: i, ii, imdp, imar, iext
    633   INTEGER :: j, jj, jmdp, jmar, nn
    634   REAL    :: xpi, zlenx, weighx, xincr,  zbordnor, zmeanor, zweinor, zbordest
    635   REAL    :: rad, zleny, weighy, masque, zbordsud, zmeasud, zweisud, zbordoue
    636 !-------------------------------------------------------------------------------
    637   imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp")
    638   jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp")
    639   imar=assert_eq(SIZE(x),SIZE(zphi,1),SIZE(mask,1),TRIM(modname)//" imar")-1
    640   jmar=assert_eq(SIZE(y),SIZE(zphi,2),SIZE(mask,2),TRIM(modname)//" jmar")
    641   IF(imar/=iim)   CALL abort_gcm(TRIM(modname),'imar/=iim'  ,1)
    642   IF(jmar/=jjm+1) CALL abort_gcm(TRIM(modname),'jmar/=jjm+1',1)
    643   iext=imdp/10
    644   xpi = ACOS(-1.)
    645   rad = 6371229.
    646 
    647 !--- ARE WE USING A READ MASK ?
    648   masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0
    649   WRITE(lunout,*)'Masque lu: ',masque_lu
    650 
    651 !--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES:
    652   ALLOCATE(xusn(imdp+2*iext))
    653   xusn(1     +iext:imdp  +iext)=xd(:)
    654   xusn(1          :       iext)=xd(1+imdp-iext:imdp)-2.*xpi
    655   xusn(1+imdp+iext:imdp+2*iext)=xd(1          :iext)+2.*xpi
    656 
    657   ALLOCATE(yusn(jmdp+2))
    658   yusn(1       )=yd(1)   +(yd(1)   -yd(2))
    659   yusn(2:jmdp+1)=yd(:)
    660   yusn(  jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1))
    661 
    662   ALLOCATE(zusn(imdp+2*iext,jmdp+2))
    663   zusn(1       +iext:imdp  +iext,2:jmdp+1)=zd  (:                   ,     :)
    664   zusn(1            :       iext,2:jmdp+1)=zd  (imdp-iext+1:imdp    ,     :)
    665   zusn(1+imdp  +iext:imdp+2*iext,2:jmdp+1)=zd  (1:iext              ,     :)
    666   zusn(1            :imdp/2+iext,       1)=zusn(1+imdp/2:imdp  +iext,     2)
    667   zusn(1+imdp/2+iext:imdp+2*iext,       1)=zusn(1       :imdp/2+iext,     2)
    668   zusn(1            :imdp/2+iext,  jmdp+2)=zusn(1+imdp/2:imdp  +iext,jmdp+1)
    669   zusn(1+imdp/2+iext:imdp+2*iext,  jmdp+2)=zusn(1       :imdp/2+iext,jmdp+1)
    670 
    671 !--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)
    672   ALLOCATE(a(imar+1),b(imar+1))
    673   b(1:imar)=(x(1:imar  )+ x(2:imar+1))/2.0
    674   b(imar+1)= x(  imar+1)+(x(  imar+1)-x(imar))/2.0
    675   a(1)=x(1)-(x(2)-x(1))/2.0
    676   a(2:imar+1)= b(1:imar)
    677 
    678   ALLOCATE(c(jmar),d(jmar))
    679   d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0
    680   d(  jmar  )= y(  jmar  )+(y(  jmar)-y(jmar-1))/2.0
    681   c(1)=y(1)-(y(2)-y(1))/2.0
    682   c(2:jmar)=d(1:jmar-1)
    683 
    684 !--- INITIALIZATIONS:
    685   ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0
    686   ALLOCATE(zmea  (imar+1,jmar)); zmea  (:,:)= 0.0
    687 
    688 !--- SUMMATION OVER GRIDPOINT AREA
    689   zleny=xpi/REAL(jmdp)*rad
    690   xincr=xpi/REAL(jmdp)/2.
    691   ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0.
    692   ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0.
    693   DO ii = 1, imar+1
    694     DO jj = 1, jmar
    695       DO j = 2,jmdp+1
    696         zlenx  =zleny  *COS(yusn(j))
    697         zbordnor=(xincr+c(jj)-yusn(j))*rad
    698         zbordsud=(xincr-d(jj)+yusn(j))*rad
    699         weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny))
    700         IF(weighy/=0) THEN
    701           DO i = 2, imdp+2*iext-1
    702             zbordest=(xusn(i)-a(ii)+xincr)*rad*COS(yusn(j))
    703             zbordoue=(b(ii)+xincr-xusn(i))*rad*COS(yusn(j))
    704             weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx))
    705             IF(weighx/=0)THEN
    706               num_tot(ii,jj)=num_tot(ii,jj)+1.0
    707               IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
    708               weight(ii,jj)=weight(ii,jj)+weighx*weighy
    709               zmea  (ii,jj)=zmea  (ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN
    710             END IF
    711           END DO
    712         END IF
    713       END DO
    714     END DO
    715   END DO
    716 
    717 !--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME
    718   IF(.NOT.masque_lu) THEN
    719     WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
    720   END IF
    721   nn=COUNT(weight(:,1:jmar-1)==0.0)
    722   IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn
    723   WHERE(weight/=0.0) zmea(:,:)=zmea(:,:)/weight(:,:)
    724 
    725 !--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD (<10%: SURF PARAMS MEANINGLESS)
    726   ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0
    727   WHERE(mask>=0.1) mask_tmp = 1.
    728   WHERE(weight(:,:)/=0.0)
    729     zphi(:,:)=mask_tmp(:,:)*zmea(:,:)
    730     zmea(:,:)=mask_tmp(:,:)*zmea(:,:)
    731   END WHERE
    732   WRITE(lunout,*)'  MEAN ORO:' ,MAXVAL(zmea)
    733 
    734 !--- Values at poles
    735   zphi(imar+1,:)=zphi(1,:)
    736 
    737   zweinor=SUM(weight(1:imar,   1),DIM=1)
    738   zweisud=SUM(weight(1:imar,jmar),DIM=1)
    739   zmeanor=SUM(weight(1:imar,   1)*zmea(1:imar,   1),DIM=1)
    740   zmeasud=SUM(weight(1:imar,jmar)*zmea(1:imar,jmar),DIM=1)
    741   zphi(:,1)=zmeanor/zweinor; zphi(:,jmar)=zmeasud/zweisud
    742 
    743 END SUBROUTINE grid_noro0
    744641!
    745642!-------------------------------------------------------------------------------

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

    r2346 r2576  
    33!*******************************************************************************
    44
     5  USE print_control_mod, ONLY: lunout
     6  USE assert_eq_m,       ONLY: assert_eq
    57  PRIVATE
    6   PUBLIC :: grid_noro
     8  PUBLIC :: grid_noro, grid_noro0
    79
    810
    … …  
    4547!    on the edge are contained in input cells.
    4648!===============================================================================
    47   USE assert_eq_m, ONLY: assert_eq
    48   USE print_control_mod, ONLY: lunout
    49   IMPLICIT NONE
    50   REAL, PARAMETER :: epsfra = 1.e-5
     49  IMPLICIT NONE
    5150!-------------------------------------------------------------------------------
    5251! Arguments:
    … …  
    315314!-------------------------------------------------------------------------------
    316315!
     316SUBROUTINE grid_noro0(xd,yd,zd,x,y,zphi,mask)
     317!
     318!===============================================================================
     319! Purpose: Extracted from grid_noro to provide geopotential height for dynamics
     320!          without any call to physics subroutines.
     321!===============================================================================
     322  IMPLICIT NONE
     323!-------------------------------------------------------------------------------
     324! 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)
     330!-------------------------------------------------------------------------------
     331! Local variables:
     332  CHARACTER(LEN=256) :: modname="grid_noro0"
     333  REAL, ALLOCATABLE :: xusn(:), yusn(:)           ! dim (imdp+2*iext) (jmdp+2)
     334  REAL, ALLOCATABLE :: zusn(:,:)                  ! dim (imdp+2*iext,jmdp+2)
     335  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)
     340  LOGICAL :: masque_lu
     341  INTEGER :: i, ii, imdp, imar, iext
     342  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
     345!-------------------------------------------------------------------------------
     346  imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp")
     347  jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp")
     348  imar=assert_eq(SIZE(x),SIZE(zphi,1),SIZE(mask,1),TRIM(modname)//" imar")-1
     349  jmar=assert_eq(SIZE(y),SIZE(zphi,2),SIZE(mask,2),TRIM(modname)//" jmar")
     350! IF(imar/=iim)   CALL abort_gcm(TRIM(modname),'imar/=iim'  ,1)
     351! IF(jmar/=jjm+1) CALL abort_gcm(TRIM(modname),'jmar/=jjm+1',1)
     352  iext=imdp/10
     353  xpi = ACOS(-1.)
     354  rad = 6371229.
     355
     356!--- ARE WE USING A READ MASK ?
     357  masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0
     358  WRITE(lunout,*)'Masque lu: ',masque_lu
     359
     360!--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES:
     361  ALLOCATE(xusn(imdp+2*iext))
     362  xusn(1     +iext:imdp  +iext)=xd(:)
     363  xusn(1          :       iext)=xd(1+imdp-iext:imdp)-2.*xpi
     364  xusn(1+imdp+iext:imdp+2*iext)=xd(1          :iext)+2.*xpi
     365
     366  ALLOCATE(yusn(jmdp+2))
     367  yusn(1       )=yd(1)   +(yd(1)   -yd(2))
     368  yusn(2:jmdp+1)=yd(:)
     369  yusn(  jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1))
     370
     371  ALLOCATE(zusn(imdp+2*iext,jmdp+2))
     372  zusn(1       +iext:imdp  +iext,2:jmdp+1)=zd  (:                   ,     :)
     373  zusn(1            :       iext,2:jmdp+1)=zd  (imdp-iext+1:imdp    ,     :)
     374  zusn(1+imdp  +iext:imdp+2*iext,2:jmdp+1)=zd  (1:iext              ,     :)
     375  zusn(1            :imdp/2+iext,       1)=zusn(1+imdp/2:imdp  +iext,     2)
     376  zusn(1+imdp/2+iext:imdp+2*iext,       1)=zusn(1       :imdp/2+iext,     2)
     377  zusn(1            :imdp/2+iext,  jmdp+2)=zusn(1+imdp/2:imdp  +iext,jmdp+1)
     378  zusn(1+imdp/2+iext:imdp+2*iext,  jmdp+2)=zusn(1       :imdp/2+iext,jmdp+1)
     379
     380!--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)
     381  ALLOCATE(a(imar+1),b(imar+1))
     382  b(1:imar)=(x(1:imar  )+ x(2:imar+1))/2.0
     383  b(imar+1)= x(  imar+1)+(x(  imar+1)-x(imar))/2.0
     384  a(1)=x(1)-(x(2)-x(1))/2.0
     385  a(2:imar+1)= b(1:imar)
     386
     387  ALLOCATE(c(jmar),d(jmar))
     388  d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0
     389  d(  jmar  )= y(  jmar  )+(y(  jmar)-y(jmar-1))/2.0
     390  c(1)=y(1)-(y(2)-y(1))/2.0
     391  c(2:jmar)=d(1:jmar-1)
     392
     393!--- INITIALIZATIONS:
     394  ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0
     395  ALLOCATE(zmea  (imar+1,jmar)); zmea  (:,:)= 0.0
     396
     397!--- SUMMATION OVER GRIDPOINT AREA
     398  zleny=xpi/REAL(jmdp)*rad
     399  xincr=xpi/REAL(jmdp)/2.
     400  ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0.
     401  ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0.
     402  DO ii = 1, imar+1
     403    DO jj = 1, jmar
     404      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
     422      END DO
     423    END DO
     424  END DO
     425
     426!--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME
     427  IF(.NOT.masque_lu) THEN
     428    WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
     429  END IF
     430  nn=COUNT(weight(:,1:jmar-1)==0.0)
     431  IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn
     432  WHERE(weight/=0.0) zmea(:,:)=zmea(:,:)/weight(:,:)
     433
     434!--- 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
     441  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
     452END SUBROUTINE grid_noro0
     453!
     454!-------------------------------------------------------------------------------
     455
     456
     457!-------------------------------------------------------------------------------
     458!
    317459SUBROUTINE MVA9(x)
    318460!
Note: See TracChangeset for help on using the changeset viewer.