Ignore:
Timestamp:
Jul 18, 2016, 9:41:10 PM (8 years ago)
Author:
Laurent Fairhead
Message:

Merged trunk changes r2545:2589 into testing branch

Location:
LMDZ5/branches/testing
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • LMDZ5/branches/testing

  • LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90

    r2471 r2594  
    1212                     prad,pg,pr,pcpp,iflag_phys)
    1313  USE dimphy, ONLY: init_dimphy
    14   USE mod_grid_phy_lmdz, ONLY: klon_glo,  & ! number of atmospheric columns (on full grid)
    15                                regular_lonlat, &  ! regular longitude-latitude grid type
    16                                nbp_lon, nbp_lat, nbp_lev
    17   USE mod_phys_lmdz_para, ONLY: klon_omp, & ! number of columns (on local omp grid)
    18                                 klon_omp_begin, & ! start index of local omp subgrid
    19                                 klon_omp_end, & ! end index of local omp subgrid
    20                                 klon_mpi_begin ! start indes of columns (on local mpi grid)
    21   USE geometry_mod, ONLY : init_geometry
     14  USE inigeomphy_mod, ONLY: inigeomphy
     15  USE mod_grid_phy_lmdz, ONLY: klon_glo ! number of atmospheric columns (on full grid)
     16  USE mod_phys_lmdz_para, ONLY: klon_omp ! number of columns (on local omp grid)
    2217  USE vertical_layers_mod, ONLY : init_vertical_layers
    2318  USE infotrac, ONLY: nqtot,nqo,nbtr,tname,ttext,type_trac,&
     
    3934  USE phystokenc_mod, ONLY: init_phystokenc
    4035  USE phyaqua_mod, ONLY: iniaqua
    41   USE physics_distribution_mod, ONLY : init_physics_distribution
    42   USE regular_lonlat_mod, ONLY : init_regular_lonlat, &
    43                                  east, west, north, south, &
    44                                  north_east, north_west, &
    45                                  south_west, south_east
    46   USE mod_interface_dyn_phys, ONLY :  init_interface_dyn_phys
    4736#ifdef INCA
    4837  USE indice_sol_mod, ONLY: nbsrf, is_oce, is_sic, is_ter, is_lic
     
    9281  CHARACTER (LEN=20) :: modname = 'iniphysiq'
    9382  CHARACTER (LEN=80) :: abort_message
    94   REAL :: total_area_phy, total_area_dyn
    9583
    96   ! boundaries, on global grid
    97   REAL,ALLOCATABLE :: boundslon_reg(:,:)
    98   REAL,ALLOCATABLE :: boundslat_reg(:,:)
    99 
    100   ! global array, on full physics grid:
    101   REAL,ALLOCATABLE :: latfi_glo(:)
    102   REAL,ALLOCATABLE :: lonfi_glo(:)
    103   REAL,ALLOCATABLE :: cufi_glo(:)
    104   REAL,ALLOCATABLE :: cvfi_glo(:)
    105   REAL,ALLOCATABLE :: airefi_glo(:)
    106   REAL,ALLOCATABLE :: boundslonfi_glo(:,:)
    107   REAL,ALLOCATABLE :: boundslatfi_glo(:,:)
    108 
    109   ! local arrays, on given MPI/OpenMP domain:
    110   REAL,ALLOCATABLE,SAVE :: latfi(:)
    111   REAL,ALLOCATABLE,SAVE :: lonfi(:)
    112   REAL,ALLOCATABLE,SAVE :: cufi(:)
    113   REAL,ALLOCATABLE,SAVE :: cvfi(:)
    114   REAL,ALLOCATABLE,SAVE :: airefi(:)
    115   REAL,ALLOCATABLE,SAVE :: boundslonfi(:,:)
    116   REAL,ALLOCATABLE,SAVE :: boundslatfi(:,:)
    117 !$OMP THREADPRIVATE (latfi,lonfi,cufi,cvfi,airefi,boundslonfi,boundslatfi)
    11884
    11985#ifndef CPP_PARA
     
    12389#endif
    12490
    125   ! Initialize Physics distibution and parameters and interface with dynamics
    126   IF (ii*jj>1) THEN ! general 3D case
    127     CALL init_physics_distribution(regular_lonlat,4, &
    128                                  nbp,ii,jj+1,nlayer,communicator)
    129   ELSE ! For 1D model
    130     CALL init_physics_distribution(regular_lonlat,4, &
    131                                  1,1,1,nlayer,communicator) 
    132   ENDIF
    133   CALL init_interface_dyn_phys
     91  ! --> initialize physics distribution, global fields and geometry
     92  ! (i.e. things in phy_common or dynphy_lonlat)
     93  CALL inigeomphy(ii,jj,nlayer, &
     94               nbp, communicator, &
     95               rlatu,rlatv, &
     96               rlonu,rlonv, &
     97               aire,cu,cv)
     98
     99  ! --> now initialize things specific to the phylmd physics package
    134100 
    135   ! init regular global longitude-latitude grid points and boundaries
    136   ALLOCATE(boundslon_reg(ii,2))
    137   ALLOCATE(boundslat_reg(jj+1,2))
    138  
    139   DO i=1,ii
    140    boundslon_reg(i,east)=rlonu(i)
    141    boundslon_reg(i,west)=rlonu(i+1)
    142   ENDDO
    143 
    144   boundslat_reg(1,north)= PI/2
    145   boundslat_reg(1,south)= rlatv(1)
    146   DO j=2,jj
    147    boundslat_reg(j,north)=rlatv(j-1)
    148    boundslat_reg(j,south)=rlatv(j)
    149   ENDDO
    150   boundslat_reg(jj+1,north)= rlatv(jj)
    151   boundslat_reg(jj+1,south)= -PI/2
    152 
    153   ! Write values in module regular_lonlat_mod
    154   CALL init_regular_lonlat(ii,jj+1, rlonv(1:ii), rlatu, &
    155                            boundslon_reg, boundslat_reg)
    156 
    157   ! Generate global arrays on full physics grid
    158   ALLOCATE(latfi_glo(klon_glo),lonfi_glo(klon_glo))
    159   ALLOCATE(cufi_glo(klon_glo),cvfi_glo(klon_glo))
    160   ALLOCATE(airefi_glo(klon_glo))
    161   ALLOCATE(boundslonfi_glo(klon_glo,4))
    162   ALLOCATE(boundslatfi_glo(klon_glo,4))
    163  
    164 
    165   IF (klon_glo>1) THEN ! general case
    166     ! North pole
    167     latfi_glo(1)=rlatu(1)
    168     lonfi_glo(1)=0.
    169     cufi_glo(1) = cu(1)
    170     cvfi_glo(1) = cv(1)
    171     boundslonfi_glo(1,north_east)=0
    172     boundslatfi_glo(1,north_east)=PI/2
    173     boundslonfi_glo(1,north_west)=2*PI
    174     boundslatfi_glo(1,north_west)=PI/2
    175     boundslonfi_glo(1,south_west)=2*PI
    176     boundslatfi_glo(1,south_west)=rlatv(1)
    177     boundslonfi_glo(1,south_east)=0
    178     boundslatfi_glo(1,south_east)=rlatv(1)
    179     DO j=2,jj
    180       DO i=1,ii
    181         k=(j-2)*ii+1+i
    182         latfi_glo(k)= rlatu(j)
    183         lonfi_glo(k)= rlonv(i)
    184         cufi_glo(k) = cu((j-1)*(ii+1)+i)
    185         cvfi_glo(k) = cv((j-1)*(ii+1)+i)
    186         boundslonfi_glo(k,north_east)=rlonu(i)
    187         boundslatfi_glo(k,north_east)=rlatv(j-1)
    188         boundslonfi_glo(k,north_west)=rlonu(i+1)
    189         boundslatfi_glo(k,north_west)=rlatv(j-1)
    190         boundslonfi_glo(k,south_west)=rlonu(i+1)
    191         boundslatfi_glo(k,south_west)=rlatv(j)
    192         boundslonfi_glo(k,south_east)=rlonu(i)
    193         boundslatfi_glo(k,south_east)=rlatv(j)
    194       ENDDO
    195     ENDDO
    196     ! South pole
    197     latfi_glo(klon_glo)= rlatu(jj+1)
    198     lonfi_glo(klon_glo)= 0.
    199     cufi_glo(klon_glo) = cu((ii+1)*jj+1)
    200     cvfi_glo(klon_glo) = cv((ii+1)*jj-ii)
    201     boundslonfi_glo(klon_glo,north_east)= 0
    202     boundslatfi_glo(klon_glo,north_east)= rlatv(jj)
    203     boundslonfi_glo(klon_glo,north_west)= 2*PI
    204     boundslatfi_glo(klon_glo,north_west)= rlatv(jj)
    205     boundslonfi_glo(klon_glo,south_west)= 2*PI
    206     boundslatfi_glo(klon_glo,south_west)= -PI/2
    207     boundslonfi_glo(klon_glo,south_east)= 0
    208     boundslatfi_glo(klon_glo,south_east)= -Pi/2
    209 
    210     ! build airefi_glo(), mesh area on physics grid
    211     CALL gr_dyn_fi(1,ii+1,jj+1,klon_glo,aire,airefi_glo)
    212     ! Poles are single points on physics grid
    213     airefi_glo(1)=sum(aire(1:ii,1))
    214     airefi_glo(klon_glo)=sum(aire(1:ii,jj+1))
    215 
    216     ! Sanity check: do total planet area match between physics and dynamics?
    217     total_area_dyn=sum(aire(1:ii,1:jj+1))
    218     total_area_phy=sum(airefi_glo(1:klon_glo))
    219     IF (total_area_dyn/=total_area_phy) THEN
    220       WRITE (lunout, *) 'iniphysiq: planet total surface discrepancy !!!'
    221       WRITE (lunout, *) '     in the dynamics total_area_dyn=', total_area_dyn
    222       WRITE (lunout, *) '  but in the physics total_area_phy=', total_area_phy
    223       IF (abs(total_area_dyn-total_area_phy)>0.00001*total_area_dyn) THEN
    224         ! stop here if the relative difference is more than 0.001%
    225         abort_message = 'planet total surface discrepancy'
    226         CALL abort_gcm(modname, abort_message, 1)
    227       ENDIF
    228     ENDIF
    229   ELSE ! klon_glo==1, running the 1D model
    230     ! just copy over input values
    231     latfi_glo(1)=rlatu(1)
    232     lonfi_glo(1)=rlonv(1)
    233     cufi_glo(1)=cu(1)
    234     cvfi_glo(1)=cv(1)
    235     airefi_glo(1)=aire(1,1)
    236     boundslonfi_glo(1,north_east)=rlonu(1)
    237     boundslatfi_glo(1,north_east)=PI/2
    238     boundslonfi_glo(1,north_west)=rlonu(2)
    239     boundslatfi_glo(1,north_west)=PI/2
    240     boundslonfi_glo(1,south_west)=rlonu(2)
    241     boundslatfi_glo(1,south_west)=rlatv(1)
    242     boundslonfi_glo(1,south_east)=rlonu(1)
    243     boundslatfi_glo(1,south_east)=rlatv(1)
    244   ENDIF ! of IF (klon_glo>1)
    245 
    246101!$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/)
    247   ! Now generate local lon/lat/cu/cv/area/bounds arrays
    248   ALLOCATE(latfi(klon_omp),lonfi(klon_omp),cufi(klon_omp),cvfi(klon_omp))
    249   ALLOCATE(airefi(klon_omp))
    250   ALLOCATE(boundslonfi(klon_omp,4))
    251   ALLOCATE(boundslatfi(klon_omp,4))
    252 
    253 
    254   offset = klon_mpi_begin - 1
    255   airefi(1:klon_omp) = airefi_glo(offset+klon_omp_begin:offset+klon_omp_end)
    256   cufi(1:klon_omp) = cufi_glo(offset+klon_omp_begin:offset+klon_omp_end)
    257   cvfi(1:klon_omp) = cvfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
    258   lonfi(1:klon_omp) = lonfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
    259   latfi(1:klon_omp) = latfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
    260   boundslonfi(1:klon_omp,:) = boundslonfi_glo(offset+klon_omp_begin:offset+klon_omp_end,:)
    261   boundslatfi(1:klon_omp,:) = boundslatfi_glo(offset+klon_omp_begin:offset+klon_omp_end,:)
    262 
    263   ! copy over local grid longitudes and latitudes
    264   CALL init_geometry(klon_omp,lonfi,latfi,boundslonfi,boundslatfi, &
    265                      airefi,cufi,cvfi)
    266102
    267103  ! copy over preff , ap(), bp(), etc
  • LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/limit_netcdf.F90

    r2542 r2594  
    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!-------------------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.