Ignore:
Timestamp:
Jul 31, 2015, 7:22:21 PM (9 years ago)
Author:
dcugnet
Message:
  • Add parallel capability for ce0l.
  • Small bug in grid_noro fixed (smoothed topography was used instead of unsmoothed one for geopotential computation at north pole).
  • Removed average of mass at poles in etat0dyn_netcdf after start_init_dyn => different results in the zoomed grid case.
  • ok_etat0=n and ok_limit=y combination now works fine (if no initial state is needed, but only limit.nc file). This required:
    • to move grid_noro0 and start_init_noro0 subroutines from etat0dyn_netcdf.F90 to limit_netcdf.F90
    • to create init_ssrf_m.F90 file, so that sub-surfaces can be initialized from limit_netcdf.F90 without any etat0*_netcdf routines call).
  • Simplified somehow the corresponding code, in particular: 1) removed obsolete flags "oldice". 2) removed flag "ibar": barycentric interpolation is used everywhere (except in start_init_subsurf, still calling grille_m - to be changed soon). 3) removed useless CPP_PHY precompilation directives, considering the possibility to run ce0l without physics is useless (ce0l is dedicated to Earth physics).
File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r2315 r2336  
    1 SUBROUTINE limit_netcdf(interbar, extrap, oldice, masque)
    2 !
    3 !-------------------------------------------------------------------------------
     1MODULE limit
     2!
     3!*******************************************************************************
    44! Author : L. Fairhead, 27/01/94
    55!-------------------------------------------------------------------------------
    … …  
    1616!  *    12/2009: D. Cugnet   (f77->f90, calendars, files from coupled runs)
    1717!-------------------------------------------------------------------------------
    18 #ifndef CPP_1D
     18
     19  USE ioipsl,             ONLY: flininfo, flinopen, flinget, flinclo,          &
     20   ioconf_calendar, ioget_calendar, lock_calendar, ioget_mon_len, ioget_year_len
     21  USE assert_eq_m,        ONLY: assert_eq
     22  USE conf_dat_m,         ONLY: conf_dat2d, conf_dat3d
     23  USE dimphy,             ONLY: klon, zmasq
     24  USE phys_state_var_mod, ONLY: pctsrf, rlon, rlat
    1925  USE control_mod
    20   USE indice_sol_mod
    21   USE dimphy
    22   USE ioipsl,             ONLY : ioget_year_len
    23   USE phys_state_var_mod, ONLY : pctsrf, rlon, rlat
    24   USE netcdf,             ONLY : NF90_OPEN,    NF90_CREATE,  NF90_CLOSE,       &
    25                    NF90_DEF_DIM, NF90_DEF_VAR, NF90_PUT_VAR, NF90_PUT_ATT,     &
    26                    NF90_NOERR,   NF90_NOWRITE, NF90_DOUBLE,  NF90_GLOBAL,      &
    27                    NF90_CLOBBER, NF90_ENDDEF,  NF90_UNLIMITED, NF90_FLOAT
    28   USE inter_barxy_m, only: inter_barxy
    29   USE netcdf95,    ONLY: nf95_def_var, nf95_put_att, nf95_put_var
    30   USE grid_atob_m, ONLY: grille_m, rugosite, sea_ice
    31   USE print_control_mod, ONLY: prt_level,lunout
    32   IMPLICIT NONE
    33 !-------------------------------------------------------------------------------
    34 ! Arguments:
    35   include "dimensions.h"
    36   include "paramet.h"
    37   LOGICAL,                    INTENT(IN) :: interbar ! barycentric interpolation
    38   LOGICAL,                    INTENT(IN) :: extrap   ! SST extrapolation flag
    39   LOGICAL,                    INTENT(IN) :: oldice   ! old way ice computation
    40   REAL, DIMENSION(iip1,jjp1), INTENT(IN) :: masque   ! land mask
    41 !-------------------------------------------------------------------------------
    42 ! Local variables:
    43   include "logic.h"
    44   include "comgeom2.h"
    45   include "comconst.h"
    46 
    47 !--- INPUT NETCDF FILES NAMES --------------------------------------------------
    48   CHARACTER(LEN=20) :: icefile, sstfile, dumstr, fnam
     26  USE init_ssrf_m,        ONLY: start_init_subsurf
     27
    4928  CHARACTER(LEN=20), PARAMETER :: &
    5029  fsst(4)=['amipbc_sst_1x1.nc   ','cpl_atm_sst.nc      ','histmth_sst.nc      '&
    … …  
    5635  vsst(4)=['tosbcs    ','SISUTESW  ','tsol_oce  ','sstk      '], &
    5736  vsic(4)=['sicbcs    ','SIICECOV  ','pourc_sic ','ci        ']
    58   CHARACTER(LEN=20), PARAMETER :: frugo='Rugos.nc            ',                &
    59                                   falbe='Albedo.nc           '
    60   CHARACTER(LEN=10), PARAMETER :: vrug='RUGOS     ', valb='ALBEDO    '
     37  CHARACTER(LEN=10), PARAMETER :: &
     38  frugo='Rugos.nc  ', falbe='Albedo.nc ', frelf='Relief.nc ',    &
     39   vrug='RUGOS     ',  valb='ALBEDO    ',  vrel='RELIEF    '
     40
     41CONTAINS
     42
     43!-------------------------------------------------------------------------------
     44!
     45SUBROUTINE limit_netcdf(masque, phis, extrap)
     46!
     47!-------------------------------------------------------------------------------
     48! Author : L. Fairhead, 27/01/94
     49!-------------------------------------------------------------------------------
     50! Purpose: Boundary conditions files building for new model using climatologies.
     51!          Both grids have to be regular.
     52!-------------------------------------------------------------------------------
     53! Note: This routine is designed to work for Earth
     54!-------------------------------------------------------------------------------
     55! Modification history:
     56!  * 23/03/1994: Z. X. Li
     57!  *    09/1999: L. Fairhead (netcdf reading in LMDZ.3.3)
     58!  *    07/2001: P. Le Van
     59!  *    11/2009: L. Guez     (ozone day & night climatos, see etat0_netcdf.F90)
     60!  *    12/2009: D. Cugnet   (f77->f90, calendars, files from coupled runs)
     61!-------------------------------------------------------------------------------
     62#ifndef CPP_1D
     63  USE indice_sol_mod
     64  USE netcdf,             ONLY: NF90_OPEN,    NF90_CREATE,  NF90_CLOSE,        &
     65                  NF90_DEF_DIM, NF90_DEF_VAR, NF90_PUT_VAR, NF90_PUT_ATT,      &
     66                  NF90_NOERR,   NF90_NOWRITE, NF90_DOUBLE,  NF90_GLOBAL,       &
     67                  NF90_CLOBBER, NF90_ENDDEF,  NF90_UNLIMITED, NF90_FLOAT
     68  USE inter_barxy_m,      ONLY: inter_barxy
     69  USE netcdf95,           ONLY: nf95_def_var, nf95_put_att, nf95_put_var
     70  IMPLICIT NONE
     71!-------------------------------------------------------------------------------
     72! Arguments:
     73  include "iniprint.h"
     74  include "dimensions.h"
     75  include "paramet.h"
     76  REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: masque ! land mask
     77  REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: phis   ! ground geopotential
     78  LOGICAL,                    INTENT(IN)    :: extrap ! SST extrapolation flag
     79!-------------------------------------------------------------------------------
     80! Local variables:
     81  include "logic.h"
     82  include "comgeom2.h"
     83  include "comconst.h"
     84
     85!--- INPUT NETCDF FILES NAMES --------------------------------------------------
     86  CHARACTER(LEN=20) :: icefile, sstfile, dumstr, fnam
    6187  CHARACTER(LEN=10) :: varname
     88
    6289!--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------
    6390  REAL               :: fi_ice(klon), verif(klon)
    … …  
    82109  CALL inigeom
    83110
     111!--- MASK, GROUND GEOPOT. & SUBSURFACES COMPUTATION (IN CASE ok_etat0==.FALSE.)
     112   IF(ALL(masque==-99999.)) THEN
     113    CALL start_init_orog0(rlonv,rlatu,phis,masque)
     114    CALL gr_dyn_fi(1,iip1,jjp1,klon,masque,zmasq)          !--- To physical grid
     115    ALLOCATE(rlon(klon),rlat(klon),pctsrf(klon,nbsrf))
     116    CALL start_init_subsurf(.FALSE.)
     117  END IF
     118
    84119!--- Beware: anneeref (from gcm.def) is used to determine output time sampling
    85120  ndays=ioget_year_len(anneeref)
    … …  
    87122!--- RUGOSITY TREATMENT --------------------------------------------------------
    88123  CALL msg(1,'Traitement de la rugosite')
    89   CALL get_2Dfield(frugo,vrug,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:))
     124  CALL get_2Dfield(frugo,vrug,'RUG',ndays,phy_rug,mask=masque(1:iim,:))
    90125
    91126!--- OCEAN TREATMENT -----------------------------------------------------------
    … …  
    108143  CALL msg(-1,'Fichier choisi pour la glace de mer:'//TRIM(icefile))
    109144
    110   CALL get_2Dfield(icefile,varname, 'SIC',interbar,ndays,phy_ice,flag=oldice)
     145  CALL get_2Dfield(icefile,varname, 'SIC',ndays,phy_ice)
    111146
    112147  ALLOCATE(pctsrf_t(klon,nbsrf,ndays))
    … …  
    167202  CALL msg(-1,'Fichier choisi pour la temperature de mer: '//TRIM(sstfile))
    168203
    169   CALL get_2Dfield(sstfile,varname,'SST',interbar,ndays,phy_sst,flag=extrap)
     204  CALL get_2Dfield(sstfile,varname,'SST',ndays,phy_sst,flag=extrap)
    170205
    171206!--- ALBEDO TREATMENT ----------------------------------------------------------
    172207  CALL msg(1,'Traitement de l albedo')
    173   CALL get_2Dfield(falbe,valb,'ALB',interbar,ndays,phy_alb)
     208  CALL get_2Dfield(falbe,valb,'ALB',ndays,phy_alb)
    174209
    175210!--- REFERENCE GROUND HEAT FLUX TREATMENT --------------------------------------
    … …  
    251286!-------------------------------------------------------------------------------
    252287!
    253 SUBROUTINE get_2Dfield(fnam, varname, mode, ibar, ndays, champo, flag, mask)
     288SUBROUTINE get_2Dfield(fnam, varname, mode, ndays, champo, flag, mask)
    254289!
    255290!-----------------------------------------------------------------------------
    … …  
    263298       NF90_CLOSE, NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION, NF90_GET_VAR, &
    264299       NF90_GET_ATT
    265   USE dimphy, ONLY : klon
    266   USE phys_state_var_mod, ONLY : pctsrf
    267   USE conf_dat_m, ONLY: conf_dat2d
    268   USE control_mod
    269300  USE pchsp_95_m, only: pchsp_95
    270301  USE pchfe_95_m, only: pchfe_95
    … …  
    281312  CHARACTER(LEN=10), INTENT(IN)     :: varname  ! NetCDF variable name
    282313  CHARACTER(LEN=3),  INTENT(IN)     :: mode     ! RUG, SIC, SST or ALB
    283   LOGICAL,           INTENT(IN)     :: ibar     ! interp on pressure levels
    284314  INTEGER,           INTENT(IN)     :: ndays    ! current year number of days
    285315  REAL,    POINTER,  DIMENSION(:, :) :: champo  ! output field = f(t)
    … …  
    312342  CHARACTER(LEN=25) :: title              ! for messages
    313343  LOGICAL           :: extrp              ! flag for extrapolation
    314   LOGICAL           :: oldice             ! flag for old way ice computation
    315344  REAL              :: chmin, chmax
    316345  INTEGER ierr
    … …  
    328357  CASE('ALB'); title='Albedo'
    329358  END SELECT
    330  
    331 
    332   extrp=.FALSE.
    333   oldice=.FALSE.
    334   IF ( PRESENT(flag) ) THEN
    335     IF ( flag .AND. mode=='SST' ) extrp=.TRUE.
    336     IF ( flag .AND. mode=='SIC' ) oldice=.TRUE.
    337   END IF
     359  extrp=.FALSE.; IF(PRESENT(flag).AND.mode=='SST') extrp=flag
    338360
    339361!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
    … …  
    400422  DO l=1, lmdep
    401423    CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,l],[imdep,jmdep,1]),fnam)
    402     CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, ibar)
     424    CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.)
    403425    IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
    404 
    405     IF(ibar .AND. .NOT.oldice) THEN
    406       IF(l==1) THEN
     426    IF(l==1) THEN
    407427      CALL msg(5,"----------------------------------------------------------")
    408428      CALL msg(5,"$$$ Interpolation barycentrique pour "//TRIM(title)//" $$$")
    409429      CALL msg(5,"----------------------------------------------------------")
    410       END IF
    411       IF(mode=='RUG') champ=LOG(champ)
    412       CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint)
    413       IF(mode=='RUG') THEN
    414         champint=EXP(champint)
    415         WHERE(NINT(mask)/=1) champint=0.001
    416       END IF
    417     ELSE
    418       SELECT CASE(mode)
    419         CASE('RUG');       CALL rugosite(dlon,dlat,champ,rlonv,rlatu,champint,mask)
    420         CASE('SIC');       CALL sea_ice (dlon,dlat,champ,rlonv,rlatu,champint)
    421         CASE('SST','ALB'); CALL grille_m(dlon,dlat,champ,rlonv,rlatu,champint)
    422       END SELECT
     430    END IF
     431    IF(mode=='RUG') champ=LOG(champ)
     432    CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint)
     433    IF(mode=='RUG') THEN
     434      champint=EXP(champint)
     435      WHERE(NINT(mask)/=1) champint=0.001
    423436    END IF
    424437    champtime(:, :, l)=champint
    … …  
    501514!-------------------------------------------------------------------------------
    502515!
     516SUBROUTINE start_init_orog0(lon_in,lat_in,phis,masque)
     517!
     518!-------------------------------------------------------------------------------
     519  IMPLICIT NONE
     520!===============================================================================
     521! Purpose:  Compute "phis" just like it would be in start_init_orog.
     522!===============================================================================
     523! Arguments:
     524  REAL,             INTENT(IN)    :: lon_in(:), lat_in(:)   ! dim (iml) (jml)
     525  REAL,             INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml)
     526!-------------------------------------------------------------------------------
     527! Local variables:
     528  CHARACTER(LEN=256) :: modname="start_init_orog0"
     529  INTEGER            :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1)
     530  REAL               :: lev(1), date, dt, deg2rad
     531  REAL, ALLOCATABLE  :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:)
     532  REAL, ALLOCATABLE  :: lat_rad(:), lat_ini(:), lat_rel(:,:)
     533!-------------------------------------------------------------------------------
     534  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml")
     535  jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml")
     536  IF(iml/=iip1) CALL abort_gcm(TRIM(modname),'iml/=iip1',1)
     537  IF(jml/=jjp1) CALL abort_gcm(TRIM(modname),'jml/=jjp1',1)
     538  pi=2.0*ASIN(1.0); deg2rad=pi/180.0
     539  IF(ANY(phis/=-99999.)) RETURN                  !--- phis ALREADY KNOWN
     540
     541!--- HIGH RESOLUTION OROGRAPHY
     542  CALL flininfo(frelf, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
     543
     544  ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
     545  CALL flinopen(frelf, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,   &
     546                lev, ttm_tmp, itau, date, dt, fid)
     547  ALLOCATE(relief_hi(iml_rel,jml_rel))
     548  CALL flinget(fid, vrel, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
     549  CALL flinclo(fid)
     550
     551!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
     552  ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
     553  lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
     554  lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
     555
     556!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
     557  ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
     558  CALL conf_dat2d(vrel, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.)
     559  DEALLOCATE(lon_ini,lat_ini)
     560
     561!--- COMPUTING SURFACE GEOPOTENTIAL USING ROUTINE grid_noro0
     562  WRITE(lunout,*)
     563  WRITE(lunout,*)'*** Compute surface geopotential ***'
     564
     565!--- CALL OROGRAPHY MODULE (REDUCED VERSION) TO COMPUTE FIELDS
     566  CALL grid_noro0(lon_rad, lat_rad, relief_hi, lon_in, lat_in, phis, masque)
     567  phis = phis * 9.81
     568  phis(iml,:) = phis(1,:)
     569  DEALLOCATE(relief_hi,lon_rad,lat_rad)
     570
     571END SUBROUTINE start_init_orog0
     572!
     573!-------------------------------------------------------------------------------
     574
     575
     576!-------------------------------------------------------------------------------
     577!
     578SUBROUTINE grid_noro0(xd,yd,zd,x,y,zphi,mask)
     579!
     580!===============================================================================
     581! Purpose: Extracted from grid_noro to provide geopotential height for dynamics
     582!          without any call to physics subroutines.
     583!===============================================================================
     584  IMPLICIT NONE
     585!-------------------------------------------------------------------------------
     586! Arguments:
     587  REAL, INTENT(IN)   :: xd(:), yd(:) !--- INPUT  COORDINATES     (imdp) (jmdp)
     588  REAL, INTENT(IN)   :: zd(:,:)      !--- INPUT  FIELD           (imdp,jmdp)
     589  REAL, INTENT(IN)   :: x(:), y(:)   !--- OUTPUT COORDINATES     (imar+1) (jmar)
     590  REAL, INTENT(OUT)  :: zphi(:,:)    !--- GEOPOTENTIAL           (imar+1,jmar)
     591  REAL, INTENT(INOUT):: mask(:,:)    !--- MASK                   (imar+1,jmar)
     592!-------------------------------------------------------------------------------
     593! Local variables:
     594  CHARACTER(LEN=256) :: modname="grid_noro0"
     595  REAL, ALLOCATABLE :: xusn(:), yusn(:)           ! dim (imdp+2*iext) (jmdp+2)
     596  REAL, ALLOCATABLE :: zusn(:,:)                  ! dim (imdp+2*iext,jmdp+2)
     597  REAL, ALLOCATABLE :: weight(:,:)                ! dim (imar+1,jmar)
     598  REAL, ALLOCATABLE :: mask_tmp(:,:), zmea(:,:)   ! dim (imar+1,jmar)
     599  REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax)
     600  REAL, ALLOCATABLE :: a(:), b(:)                 ! dim (imax)
     601  REAL, ALLOCATABLE :: c(:), d(:)                 ! dim (jmax)
     602  LOGICAL :: masque_lu
     603  INTEGER :: i, ii, imdp, imar, iext
     604  INTEGER :: j, jj, jmdp, jmar, nn
     605  REAL    :: xpi, zlenx, weighx, xincr,  zbordnor, zmeanor, zweinor, zbordest
     606  REAL    :: rad, zleny, weighy, masque, zbordsud, zmeasud, zweisud, zbordoue
     607!-------------------------------------------------------------------------------
     608  imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp")
     609  jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp")
     610  imar=assert_eq(SIZE(x),SIZE(zphi,1),SIZE(mask,1),TRIM(modname)//" imar")-1
     611  jmar=assert_eq(SIZE(y),SIZE(zphi,2),SIZE(mask,2),TRIM(modname)//" jmar")
     612  IF(imar/=iim)   CALL abort_gcm(TRIM(modname),'imar/=iim'  ,1)
     613  IF(jmar/=jjm+1) CALL abort_gcm(TRIM(modname),'jmar/=jjm+1',1)
     614  iext=imdp/10
     615  xpi = ACOS(-1.)
     616  rad = 6371229.
     617
     618!--- ARE WE USING A READ MASK ?
     619  masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0
     620  WRITE(lunout,*)'Masque lu: ',masque_lu
     621
     622!--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES:
     623  ALLOCATE(xusn(imdp+2*iext))
     624  xusn(1     +iext:imdp  +iext)=xd(:)
     625  xusn(1          :       iext)=xd(1+imdp-iext:imdp)-2.*xpi
     626  xusn(1+imdp+iext:imdp+2*iext)=xd(1          :iext)+2.*xpi
     627
     628  ALLOCATE(yusn(jmdp+2))
     629  yusn(1       )=yd(1)   +(yd(1)   -yd(2))
     630  yusn(2:jmdp+1)=yd(:)
     631  yusn(  jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1))
     632
     633  ALLOCATE(zusn(imdp+2*iext,jmdp+2))
     634  zusn(1       +iext:imdp  +iext,2:jmdp+1)=zd  (:                   ,     :)
     635  zusn(1            :       iext,2:jmdp+1)=zd  (imdp-iext+1:imdp    ,     :)
     636  zusn(1+imdp  +iext:imdp+2*iext,2:jmdp+1)=zd  (1:iext              ,     :)
     637  zusn(1            :imdp/2+iext,       1)=zusn(1+imdp/2:imdp  +iext,     2)
     638  zusn(1+imdp/2+iext:imdp+2*iext,       1)=zusn(1       :imdp/2+iext,     2)
     639  zusn(1            :imdp/2+iext,  jmdp+2)=zusn(1+imdp/2:imdp  +iext,jmdp+1)
     640  zusn(1+imdp/2+iext:imdp+2*iext,  jmdp+2)=zusn(1       :imdp/2+iext,jmdp+1)
     641
     642!--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)
     643  ALLOCATE(a(imar+1),b(imar+1))
     644  b(1:imar)=(x(1:imar  )+ x(2:imar+1))/2.0
     645  b(imar+1)= x(  imar+1)+(x(  imar+1)-x(imar))/2.0
     646  a(1)=x(1)-(x(2)-x(1))/2.0
     647  a(2:imar+1)= b(1:imar)
     648
     649  ALLOCATE(c(jmar),d(jmar))
     650  d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0
     651  d(  jmar  )= y(  jmar  )+(y(  jmar)-y(jmar-1))/2.0
     652  c(1)=y(1)-(y(2)-y(1))/2.0
     653  c(2:jmar)=d(1:jmar-1)
     654
     655!--- INITIALIZATIONS:
     656  ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0
     657  ALLOCATE(zmea  (imar+1,jmar)); zmea  (:,:)= 0.0
     658
     659!--- SUMMATION OVER GRIDPOINT AREA
     660  zleny=xpi/REAL(jmdp)*rad
     661  xincr=xpi/REAL(jmdp)/2.
     662  ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0.
     663  ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0.
     664  DO ii = 1, imar+1
     665    DO jj = 1, jmar
     666      DO j = 2,jmdp+1
     667        zlenx  =zleny  *COS(yusn(j))
     668        zbordnor=(xincr+c(jj)-yusn(j))*rad
     669        zbordsud=(xincr-d(jj)+yusn(j))*rad
     670        weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny))
     671        IF(weighy/=0) THEN
     672          DO i = 2, imdp+2*iext-1
     673            zbordest=(xusn(i)-a(ii)+xincr)*rad*COS(yusn(j))
     674            zbordoue=(b(ii)+xincr-xusn(i))*rad*COS(yusn(j))
     675            weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx))
     676            IF(weighx/=0)THEN
     677              num_tot(ii,jj)=num_tot(ii,jj)+1.0
     678              IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
     679              weight(ii,jj)=weight(ii,jj)+weighx*weighy
     680              zmea  (ii,jj)=zmea  (ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN
     681            END IF
     682          END DO
     683        END IF
     684      END DO
     685    END DO
     686  END DO
     687
     688!--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME
     689  IF(.NOT.masque_lu) THEN
     690    WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
     691  END IF
     692  nn=COUNT(weight(:,1:jmar-1)==0.0)
     693  IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn
     694  WHERE(weight/=0.0) zmea(:,:)=zmea(:,:)/weight(:,:)
     695
     696!--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD (<10%: SURF PARAMS MEANINGLESS)
     697  ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0
     698  WHERE(mask>=0.1) mask_tmp = 1.
     699  WHERE(weight(:,:)/=0.0)
     700    zphi(:,:)=mask_tmp(:,:)*zmea(:,:)
     701    zmea(:,:)=mask_tmp(:,:)*zmea(:,:)
     702  END WHERE
     703  WRITE(lunout,*)'  MEAN ORO:' ,MAXVAL(zmea)
     704
     705!--- Values at poles
     706  zphi(imar+1,:)=zphi(1,:)
     707
     708  zweinor=SUM(weight(1:imar,   1),DIM=1)
     709  zweisud=SUM(weight(1:imar,jmar),DIM=1)
     710  zmeanor=SUM(weight(1:imar,   1)*zmea(1:imar,   1),DIM=1)
     711  zmeasud=SUM(weight(1:imar,jmar)*zmea(1:imar,jmar),DIM=1)
     712  zphi(:,1)=zmeanor/zweinor; zphi(:,jmar)=zmeasud/zweisud
     713
     714END SUBROUTINE grid_noro0
     715!
     716!-------------------------------------------------------------------------------
     717
     718
     719!-------------------------------------------------------------------------------
     720!
    503721FUNCTION year_len(y,cal_in)
    504722!
    505723!-------------------------------------------------------------------------------
    506   USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_year_len
    507724  IMPLICIT NONE
    508725!-------------------------------------------------------------------------------
    … …  
    537754!
    538755!-------------------------------------------------------------------------------
    539   USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_mon_len
    540756  IMPLICIT NONE
    541757!-------------------------------------------------------------------------------
    … …  
    624840!-------------------------------------------------------------------------------
    625841  USE netcdf, ONLY : NF90_NOERR, NF90_STRERROR
    626   USE print_control_mod, ONLY: lunout
    627842  IMPLICIT NONE
    628843!-------------------------------------------------------------------------------
    … …  
    643858! of #ifndef CPP_1D
    644859END SUBROUTINE limit_netcdf
     860
     861END MODULE limit
     862!
     863!*******************************************************************************
Note: See TracChangeset for help on using the changeset viewer.