Ignore:
Timestamp:
Feb 23, 2010, 10:29:54 PM (14 years ago)
Author:
Laurent Fairhead
Message:
  • Modifications to the start and limit creation routines to account for different

calendars

  • Modification to phyetat0 to force the mask read in the start files to match the

surface fractions read in the limit file

  • Force readaerosol.F90 to read in aerosols file with 12 timesteps
  • Modifications aux routines de créations des fichiers start et limit pour prendre

en compte différents calendriers

  • Modification à phyetat0 pour forcer le masque lu dans le fichier start à être

compatible avec les fractions de surface lu dans le fichier limit

  • Forcer readaerosol à ne lire que des fichiers à 12 pas de temps
File:
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ4/trunk/libf/dyn3d/limit_netcdf.F90

    r1318 r1319  
    11!
    22! $Id$
    3 !
    4 C
    5 C
    6       SUBROUTINE limit_netcdf(interbar, extrap, oldice, masque)
     3!-------------------------------------------------------------------------------
     4!
     5SUBROUTINE limit_netcdf(interbar, extrap, oldice, masque)
     6!
     7!-------------------------------------------------------------------------------
     8! Author : L. Fairhead, 27/01/94
     9!-------------------------------------------------------------------------------
     10! Purpose: Boundary conditions files building for new model using climatologies.
     11!          Both grids have to be regular.
     12!-------------------------------------------------------------------------------
     13! Note: This routine is designed to work for Earth
     14!-------------------------------------------------------------------------------
     15! Modification history:
     16!  * 23/03/1994: Z. X. Li
     17!  *    09/1999: L. Fairhead (netcdf reading in LMDZ.3.3)
     18!  *    07/2001: P. Le Van
     19!  *    11/2009: L. Guez     (ozone day & night climatos, see etat0_netcdf.F90)
     20!  *    12/2009: D. Cugnet   (f77->f90, calendars, files from coupled runs)
     21!-------------------------------------------------------------------------------
    722#ifdef CPP_EARTH
    8 ! This routine is designed to work for Earth
    9       USE dimphy
    10       use phys_state_var_mod , ONLY : pctsrf
    11       IMPLICIT none
    12 c
    13 c-------------------------------------------------------------
    14 C Author : L. Fairhead
    15 C Date   : 27/01/94
    16 C Objet  : Construction des fichiers de conditions aux limites
    17 C          pour le nouveau
    18 C          modele a partir de fichiers de climatologie. Les deux
    19 C          grilles doivent etre regulieres
    20 c
    21 c Modifie par z.x.li (le23mars1994)
    22 c Modifie par L. Fairhead (fairhead@lmd.jussieu.fr) septembre 1999
    23 c                         pour lecture netcdf dans LMDZ.3.3
    24 c Modifie par P;Le Van  ,  juillet 2001
    25 c-------------------------------------------------------------
    26 c
     23  USE dimphy
     24  USE ioipsl,             ONLY : ioget_year_len
     25  USE phys_state_var_mod, ONLY : pctsrf
     26  USE netcdf,             ONLY : NF90_OPEN,    NF90_CREATE,  NF90_CLOSE,       &
     27                   NF90_DEF_DIM, NF90_DEF_VAR, NF90_PUT_VAR, NF90_PUT_ATT,     &
     28                   NF90_NOERR,   NF90_NOWRITE, NF90_DOUBLE,  NF90_GLOBAL,      &
     29                   NF90_CLOBBER, NF90_ENDDEF,  NF90_UNLIMITED
     30#endif
     31  IMPLICIT NONE
     32!-------------------------------------------------------------------------------
     33! Arguments:
    2734#include "dimensions.h"
    2835#include "paramet.h"
     36#include "iniprint.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#ifndef CPP_EARTH
     42  WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
     43#else
     44!-------------------------------------------------------------------------------
     45! Local variables:
    2946#include "control.h"
    3047#include "logic.h"
    31 #include "netcdf.inc"
    3248#include "comvert.h"
    3349#include "comgeom2.h"
    3450#include "comconst.h"
    35 cy#include "dimphy.h"
     51#include "indicesol.h"
     52
     53!--- For fractionary sub-cell use (old coding used soil index: 0,1,2,3) --------
     54  LOGICAL, PARAMETER :: fracterre=.TRUE.
     55
     56!--- INPUT NETCDF FILES NAMES --------------------------------------------------
     57  CHARACTER(LEN=25) :: icefile, sstfile, dumstr
     58  CHARACTER(LEN=25), PARAMETER :: famipsst='amipbc_sst_1x1.nc        ',        &
     59                                  famipsic='amipbc_sic_1x1.nc        ',        &
     60                                  fclimsst='amipbc_sst_1x1_clim.nc   ',        &
     61                                  fclimsic='amipbc_sic_1x1_clim.nc   ',        &
     62                                  fcpldsst='cpl_atm_sst.nc           ',        &
     63                                  fcpldsic='cpl_atm_sic.nc           ',        &
     64                                  frugo   ='Rugos.nc                 ',        &
     65                                  falbe   ='Albedo.nc                '
     66
     67!--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------
     68  REAL,   DIMENSION(klon)                :: fi_ice, verif
     69  REAL,   DIMENSION(:,:),   POINTER      :: phy_rug=>NULL(), phy_ice=>NULL()
     70  REAL,   DIMENSION(:,:),   POINTER      :: phy_sst=>NULL(), phy_alb=>NULL()
     71  REAL,   DIMENSION(:,:),   ALLOCATABLE  :: phy_bil
     72  REAL,   DIMENSION(:,:,:), ALLOCATABLE  :: pctsrf_t
     73  INTEGER                                :: nbad
     74
     75!--- VARIABLES FOR OUTPUT FILE WRITING -----------------------------------------
     76  INTEGER :: ierr, nid, ndim, ntim, k
     77  INTEGER, DIMENSION(2) :: dims
     78  INTEGER :: id_tim,  id_SST,  id_BILS, id_RUG, id_ALB
     79  INTEGER :: id_FOCE, id_FSIC, id_FTER, id_FLIC
     80  INTEGER :: NF90_FORMAT
     81  LOGICAL :: lCPL                    !--- T: IPCC-IPSL cpl model output files
     82
     83!--- MICS (PHYSICAL KEYS, TIME) ------------------------------------------------
     84!  INTEGER, PARAMETER        :: longcles=20
     85!  REAL, DIMENSION(longcles) :: clesphy0
     86  INTEGER                   :: ndays
     87 
     88!--- INITIALIZATIONS -----------------------------------------------------------
     89#ifdef NC_DOUBLE
     90  NF90_FORMAT=NF90_DOUBLE
     91#else
     92  NF90_FORMAT=NF90_FLOAT
     93#endif
     94
     95!  CALL conf_gcm(99, .TRUE., clesphy0)
     96  pi    = 4.*ATAN(1.)
     97  rad   = 6371229.
     98  daysec= 86400.
     99  omeg  = 2.*pi/daysec
     100  g     = 9.8
     101  kappa = 0.2857143
     102  cpp   = 1004.70885
     103  dtvr  = daysec/FLOAT(day_step)
     104  CALL inigeom
     105
     106!--- Beware: anneeref (from gcm.def) is used to determine output time sampling
     107  ndays=ioget_year_len(anneeref)
     108
     109!--- RUGOSITY TREATMENT --------------------------------------------------------
     110  WRITE(lunout,*) 'Traitement de la rugosite'
     111  CALL get_2Dfield(frugo,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:))
     112
     113!--- OCEAN TREATMENT -----------------------------------------------------------
     114  PRINT*, 'Traitement de la glace oceanique' ; icefile=''; lCPL=.FALSE.
     115
     116! Input SIC file selection
     117  icefile='fake'
     118  IF(NF90_OPEN(famipsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(famipsic)
     119  IF(NF90_OPEN(fclimsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fclimsic)
     120  IF(NF90_OPEN(fcpldsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fcpldsic)
     121  SELECT CASE(icefile)
     122    CASE(famipsic); dumstr='Amip.'
     123    CASE(fclimsic); dumstr='Amip climatologique.'
     124    CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE.
     125    CASE('fake');   CALL abort_gcm('limit_netcdf','Fichier SIC non reconnu.',1)
     126  END SELECT
     127  ierr=NF90_CLOSE(nid)
     128  WRITE(lunout,*)'Pour la glace de mer a ete choisi un fichier '//TRIM(dumstr)
     129
     130  CALL get_2Dfield(icefile,'SIC',interbar,ndays,phy_ice,flag=oldice)
     131
     132  ALLOCATE(pctsrf_t(klon,nbsrf,ndays))
     133  DO k=1,ndays
     134    fi_ice=phy_ice(:,k)
     135    WHERE(fi_ice>=1.0  ) fi_ice=1.0
     136    WHERE(fi_ice<EPSFRA) fi_ice=0.0
     137    IF(fracterre) THEN
     138      pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)       ! land soil
     139      pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic)       ! land ice
     140      IF(lCPL) THEN                              ! SIC=pICE*(1-LIC-TER)
     141        pctsrf_t(:,is_sic,k)=fi_ice*(1-pctsrf(:,is_lic)-pctsrf(:,is_ter))
     142      ELSE                                        ! SIC=pICE-LIC
     143        pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k)
     144      END IF
     145      WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0.
     146      WHERE(1.0-zmasq<EPSFRA)
     147        pctsrf_t(:,is_sic,k)=0.0
     148        pctsrf_t(:,is_oce,k)=0.0
     149      ELSEWHERE
     150        WHERE(pctsrf_t(:,is_sic,k)>=1.0-zmasq)
     151          pctsrf_t(:,is_sic,k)=1.0-zmasq
     152          pctsrf_t(:,is_oce,k)=0.0
     153        ELSEWHERE
     154          pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)
     155          WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)
     156            pctsrf_t(:,is_oce,k)=0.0
     157            pctsrf_t(:,is_sic,k)=1.0-zmasq
     158          END WHERE
     159        END WHERE
     160      END WHERE
     161      nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0)
     162      IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
     163      nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA)
     164      IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad
     165    ELSE
     166      pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)
     167      WHERE(NINT(pctsrf(:,is_ter))==1)
     168        pctsrf_t(:,is_sic,k)=0.
     169        pctsrf_t(:,is_oce,k)=0.                 
     170        WHERE(fi_ice>=1.e-5)
     171          pctsrf_t(:,is_lic,k)=fi_ice
     172          pctsrf_t(:,is_ter,k)=pctsrf_t(:,is_ter,k)-pctsrf_t(:,is_lic,k)
     173        ELSEWHERE
     174          pctsrf_t(:,is_lic,k)=0.0
     175        END WHERE
     176      ELSEWHERE
     177        pctsrf_t(:,is_lic,k) = 0.0
     178        WHERE(fi_ice>=1.e-5)
     179          pctsrf_t(:,is_ter,k)=0.0
     180          pctsrf_t(:,is_sic,k)=fi_ice
     181          pctsrf_t(:,is_oce,k)=1.0-pctsrf_t(:,is_sic,k)
     182        ELSEWHERE
     183          pctsrf_t(:,is_sic,k)=0.0
     184          pctsrf_t(:,is_oce,k)=1.0
     185        END WHERE
     186      END WHERE
     187      verif=sum(pctsrf_t(:,:,k),dim=2)
     188      nbad=COUNT(verif<1.0-1.e-5.OR.verif>1.0+1.e-5)
     189      IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
     190    END IF
     191  END DO
     192  DEALLOCATE(phy_ice)
     193
     194!--- SST TREATMENT -------------------------------------------------------------
     195  WRITE(lunout,*) 'Traitement de la sst' ; sstfile=''; lCPL=.FALSE.
     196
     197! Input SST file selection
     198  sstfile='fake'
     199  IF(NF90_OPEN(famipsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(famipsst)
     200  IF(NF90_OPEN(fclimsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fclimsst)
     201  IF(NF90_OPEN(fcpldsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fcpldsst)
     202  SELECT CASE(icefile)
     203    CASE(famipsic); dumstr='Amip.'
     204    CASE(fclimsic); dumstr='Amip climatologique.'
     205    CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE.
     206    CASE('fake');   CALL abort_gcm('limit_netcdf','Fichier SST non reconnu',1)
     207  END SELECT
     208  ierr=NF90_CLOSE(nid)
     209  WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(dumstr)
     210
     211  CALL get_2Dfield(trim(sstfile),'SST',interbar,ndays,phy_sst,flag=extrap)
     212
     213!--- ALBEDO TREATMENT ----------------------------------------------------------
     214  WRITE(lunout,*) 'Traitement de l albedo'
     215  CALL get_2Dfield(falbe,'ALB',interbar,ndays,phy_alb)
     216
     217!--- REFERENCE GROUND HEAT FLUX TREATMENT --------------------------------------
     218  ALLOCATE(phy_bil(klon,ndays)); phy_bil=0.0
     219
     220!--- OUTPUT FILE WRITING -------------------------------------------------------
     221  WRITE(lunout,*) 'Ecriture du fichier limit'
     222
     223  !--- File creation
     224  ierr=NF90_CREATE("limit.nc",NF90_CLOBBER,nid)
     225  ierr=NF90_PUT_ATT(nid,NF90_GLOBAL,"title","Fichier conditions aux limites")
     226
     227  !--- Dimensions creation
     228  ierr=NF90_DEF_DIM(nid,"points_physiques",klon,ndim)
     229  ierr=NF90_DEF_DIM(nid,"time",NF90_UNLIMITED,ntim)
     230
     231  dims=(/ndim,ntim/)
     232
     233  !--- Variables creation
     234  ierr=NF90_DEF_VAR(nid,"TEMPS",NF90_FORMAT,dims,id_tim)
     235  ierr=NF90_DEF_VAR(nid,"FOCE", NF90_FORMAT,dims,id_FOCE)
     236  ierr=NF90_DEF_VAR(nid,"FSIC", NF90_FORMAT,dims,id_FSIC)
     237  ierr=NF90_DEF_VAR(nid,"FTER", NF90_FORMAT,dims,id_FTER)
     238  ierr=NF90_DEF_VAR(nid,"FLIC", NF90_FORMAT,dims,id_FLIC)
     239  ierr=NF90_DEF_VAR(nid,"SST",  NF90_FORMAT,dims,id_SST)
     240  ierr=NF90_DEF_VAR(nid,"BILS", NF90_FORMAT,dims,id_BILS)
     241  ierr=NF90_DEF_VAR(nid,"ALB",  NF90_FORMAT,dims,id_ALB)
     242  ierr=NF90_DEF_VAR(nid,"RUG",  NF90_FORMAT,dims,id_RUG)
     243
     244  !--- Attributes creation
     245  ierr=NF90_PUT_ATT(nid,id_tim, "title","Jour dans l annee")
     246  ierr=NF90_PUT_ATT(nid,id_FOCE,"title","Fraction ocean")
     247  ierr=NF90_PUT_ATT(nid,id_FSIC,"title","Fraction glace de mer")
     248  ierr=NF90_PUT_ATT(nid,id_FTER,"title","Fraction terre")
     249  ierr=NF90_PUT_ATT(nid,id_FLIC,"title","Fraction land ice")
     250  ierr=NF90_PUT_ATT(nid,id_SST ,"title","Temperature superficielle de la mer")
     251  ierr=NF90_PUT_ATT(nid,id_BILS,"title","Reference flux de chaleur au sol")
     252  ierr=NF90_PUT_ATT(nid,id_ALB, "title","Albedo a la surface")
     253  ierr=NF90_PUT_ATT(nid,id_RUG, "title","Rugosite")
     254
     255  ierr=NF90_ENDDEF(nid)
     256
     257  !--- Variables saving
     258  ierr=NF90_PUT_VAR(nid,id_tim,(/(DBLE(k),k=1,ndays)/))
     259  ierr=NF90_PUT_VAR(nid,id_FOCE,pctsrf_t(:,is_oce,:),(/1,1/),(/klon,ndays/))
     260  ierr=NF90_PUT_VAR(nid,id_FSIC,pctsrf_t(:,is_sic,:),(/1,1/),(/klon,ndays/))
     261  ierr=NF90_PUT_VAR(nid,id_FTER,pctsrf_t(:,is_ter,:),(/1,1/),(/klon,ndays/))
     262  ierr=NF90_PUT_VAR(nid,id_FLIC,pctsrf_t(:,is_lic,:),(/1,1/),(/klon,ndays/))
     263  ierr=NF90_PUT_VAR(nid,id_SST ,phy_sst(:,:),(/1,1/),(/klon,ndays/))
     264  ierr=NF90_PUT_VAR(nid,id_BILS,phy_bil(:,:),(/1,1/),(/klon,ndays/))
     265  ierr=NF90_PUT_VAR(nid,id_ALB ,phy_alb(:,:),(/1,1/),(/klon,ndays/))
     266  ierr=NF90_PUT_VAR(nid,id_RUG ,phy_rug(:,:),(/1,1/),(/klon,ndays/))
     267
     268  ierr=NF90_CLOSE(nid)
     269
     270  DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug)
     271#endif
     272! of #ifdef CPP_EARTH
     273  STOP
     274
     275
     276!===============================================================================
     277!
     278  CONTAINS
     279!
     280!===============================================================================
     281
     282
     283!-------------------------------------------------------------------------------
     284!
     285SUBROUTINE get_2Dfield(fnam, mode, ibar, ndays, champo, flag, mask)
     286!
     287!-------------------------------------------------------------------------------
     288! Comments:
     289!   There are two assumptions concerning the NetCDF files, that are satisfied
     290!   with files that are conforming NC convention:
     291!     1) The last dimension of the variables used is the time record.
     292!     2) Dimensional variables have the same names as corresponding dimensions.
     293!-------------------------------------------------------------------------------
     294  USE netcdf, ONLY : NF90_OPEN,    NF90_INQ_VARID, NF90_INQUIRE_VARIABLE,      &
     295                     NF90_CLOSE,   NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION,     &
     296                     NF90_GET_VAR, NF90_GET_ATT
     297  USE dimphy, ONLY : klon
     298  USE phys_state_var_mod, ONLY : pctsrf
     299  IMPLICIT NONE
     300#include "dimensions.h"
     301#include "paramet.h"
     302#include "comgeom2.h"
     303#include "control.h"
    36304#include "indicesol.h"
    37305#include "iniprint.h"
    38 c
    39 c-----------------------------------------------------------------------
    40       LOGICAL interbar, extrap, oldice
    41 
    42       REAL phy_nat(klon,360), phy_nat0(klon)
    43       REAL phy_alb(klon,360)
    44       REAL phy_sst(klon,360)
    45       REAL phy_bil(klon,360)
    46       REAL phy_rug(klon,360)
    47       REAL phy_ice(klon)
    48 c
    49       real pctsrf_t(klon,nbsrf,360)
    50 
    51       REAL verif
    52 
    53       REAL masque(iip1,jjp1)
    54       REAL mask(iim,jjp1)
    55 CPB
    56 C newlmt indique l'utilisation de la sous-maille fractionnelle
    57 C tandis que l'ancien codage utilise l'indicateur du sol (0,1,2,3)
    58       LOGICAL newlmt, fracterre
    59       PARAMETER(newlmt=.TRUE.)
    60       PARAMETER(fracterre = .TRUE.)
    61 
    62 C Declarations pour le champ de depart
    63       INTEGER imdep, jmdep,lmdep
    64       INTEGER  tbid
    65       PARAMETER ( tbid = 60 )        ! >52 semaines
    66       REAL  timecoord(tbid)
    67 c
    68       REAL , ALLOCATABLE :: dlon_msk(:), dlat_msk(:)
    69       REAL , ALLOCATABLE :: lonmsk_ini(:), latmsk_ini(:)
    70       REAL , ALLOCATABLE :: dlon(:), dlat(:)
    71       REAL , ALLOCATABLE :: dlon_ini(:), dlat_ini(:)
    72       REAL , ALLOCATABLE :: champ_msk(:), champ(:)
    73       REAL , ALLOCATABLE :: work(:,:)
    74 
    75       CHARACTER*25 title
    76 
    77 C Declarations pour le champ interpole 2D
    78       REAL champint(iim,jjp1)
    79       real chmin,chmax
    80 
    81 C Declarations pour le champ interpole 3D
    82       REAL champtime(iim,jjp1,tbid)
    83       REAL timeyear(tbid)
    84       REAL champan(iip1,jjp1,366)
    85 
    86 C Declarations pour l'inteprolation verticale
    87       REAL ax(tbid), ay(tbid)
    88       REAL by
    89       REAL yder(tbid)
    90 
    91 
    92       INTEGER ierr
    93       INTEGER dimfirst(3)
    94       INTEGER dimlast(3)
    95 c
    96       INTEGER nid, ndim, ntim
    97       INTEGER dims(2), debut(2), epais(2)
    98       INTEGER id_tim
    99       INTEGER id_NAT, id_SST, id_BILS, id_RUG, id_ALB
    100 CPB
    101       INTEGER id_FOCE, id_FSIC, id_FTER, id_FLIC
    102 
    103       INTEGER i, j, k, l, ji
    104 c declarations pour lecture glace de mer
    105       INTEGER :: iml_lic, jml_lic, llm_tmp, ttm_tmp, iret
    106       INTEGER :: itaul(1), fid
    107       REAL :: lev(1), date, dt
    108       REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_lic, lat_lic
    109       REAL, ALLOCATABLE, DIMENSION(:)  :: dlon_lic, dlat_lic
    110       REAL, ALLOCATABLE, DIMENSION (:,:) :: fraclic
    111       REAL :: flic_tmp(iip1, jjp1)
    112 
    113 c Diverses variables locales
    114       REAL time
    115 ! pour la lecture du fichier masque ocean
    116       integer :: nid_o2a
    117       logical :: couple = .false.
    118       INTEGER :: iml_omask, jml_omask
    119       REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_omask, lat_omask
    120       REAL, ALLOCATABLE, DIMENSION(:)  :: dlon_omask, dlat_omask
    121       REAL, ALLOCATABLE, DIMENSION (:,:) :: ocemask, ocetmp
    122       real, dimension(klon) :: ocemask_fi
    123 
    124       INTEGER          longcles
    125       PARAMETER      ( longcles = 20 )
    126       REAL  clesphy0 ( longcles      )
    127 #include "serre.h"
    128       INTEGER ncid,varid,ndimid(4),dimid
    129       character*30 namedim
    130       CHARACTER*80 :: varname
    131 
    132 cIM28/02/2002 <== PM
    133       REAL tmidmonth(12)
    134       SAVE tmidmonth
    135       DATA tmidmonth/15,45,75,105,135,165,195,225,255,285,315,345/
    136 
    137 c initialisations:
    138       CALL conf_gcm( 99, .TRUE. , clesphy0 )
    139 
    140 
    141       pi     = 4. * ATAN(1.)
    142       rad    = 6 371 229.
    143       omeg   = 4.* ASIN(1.)/(24.*3600.)
    144       g      = 9.8
    145       daysec = 86400.
    146       kappa  = 0.2857143
    147       cpp    = 1004.70885
    148       dtvr    = daysec/FLOAT(day_step)
    149       CALL inigeom
    150 c
    151 C Traitement du relief au sol
    152 c
    153       write(*,*) 'Traitement du relief au sol pour fabriquer masque'
    154       ierr = NF_OPEN('Relief.nc', NF_NOWRITE, ncid)
    155 
    156       if (ierr.ne.0) then
    157         print *, NF_STRERROR(ierr)
    158         STOP
    159       ENDIF
    160 
    161       ierr = NF_INQ_VARID(ncid,'RELIEF',varid)
    162       if (ierr.ne.0) then
    163         print *, NF_STRERROR(ierr)
    164         STOP
    165       ENDIF
    166       ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
    167       if (ierr.ne.0) then
    168         print *, NF_STRERROR(ierr)
    169         STOP
    170       ENDIF
    171       ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
    172       if (ierr.ne.0) then
    173         print *, NF_STRERROR(ierr)
    174         STOP
    175       ENDIF
    176       print*,'variable ', namedim, 'dimension ', imdep
    177       ierr = NF_INQ_VARID(ncid,namedim,dimid)
    178       if (ierr.ne.0) then
    179         print *, NF_STRERROR(ierr)
    180         STOP
    181       ENDIF
    182 
    183       ALLOCATE( lonmsk_ini(imdep) )
    184       ALLOCATE(   dlon_msk(imdep) )
    185 
    186 #ifdef NC_DOUBLE
    187       ierr = NF_GET_VAR_DOUBLE(ncid,dimid,lonmsk_ini)
    188 #else
    189       ierr = NF_GET_VAR_REAL(ncid,dimid,lonmsk_ini)
    190 #endif
    191 
    192 c
    193       if (ierr.ne.0) then
    194         print *, NF_STRERROR(ierr)
    195         STOP
    196       ENDIF
    197       ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
    198       if (ierr.ne.0) then
    199         print *, NF_STRERROR(ierr)
    200         STOP
    201       ENDIF
    202       print*,'variable ', namedim, 'dimension ', jmdep
    203       ierr = NF_INQ_VARID(ncid,namedim,dimid)
    204       if (ierr.ne.0) then
    205         print *, NF_STRERROR(ierr)
    206         STOP
    207       ENDIF
    208 
    209       ALLOCATE( latmsk_ini(jmdep) )
    210       ALLOCATE(   dlat_msk(jmdep) )
    211       ALLOCATE(  champ_msk(imdep*jmdep) )
    212 
    213 #ifdef NC_DOUBLE
    214       ierr = NF_GET_VAR_DOUBLE(ncid,dimid,latmsk_ini)
    215 #else
    216       ierr = NF_GET_VAR_REAL(ncid,dimid,latmsk_ini)
    217 #endif
    218 c
    219       if (ierr.ne.0) then
    220         print *, NF_STRERROR(ierr)
    221         STOP
    222       ENDIF
    223 #ifdef NC_DOUBLE
    224       ierr = NF_GET_VAR_DOUBLE(ncid,varid,champ_msk)
    225 #else
    226       ierr = NF_GET_VAR_REAL(ncid,varid,champ_msk)
    227 #endif
    228 c
    229       if (ierr.ne.0) then
    230         print *, NF_STRERROR(ierr)
    231         STOP
    232       ENDIF
    233 c
    234       title='RELIEF'
    235 
    236       CALL conf_dat2d(title,imdep, jmdep, lonmsk_ini, latmsk_ini,
    237      . dlon_msk, dlat_msk, champ_msk, interbar  )
    238 
    239       DO i = 1, iim
    240       DO j = 1, jjp1
    241          mask(i,j) = masque(i,j)
    242       ENDDO
    243       ENDDO
    244       WRITE(*,*) 'MASK:'
    245       WRITE(*,'(96i1)')INT(mask)     
    246       ierr = NF_CLOSE(ncid)
    247 c
    248 c
    249 C Traitement de la rugosite
    250 c
    251       PRINT*, 'Traitement de la rugosite'
    252       ierr = NF_OPEN('Rugos.nc', NF_NOWRITE, ncid)
    253       if (ierr.ne.0) then
    254         print *, NF_STRERROR(ierr)
    255         STOP
    256       ENDIF
    257 
    258       ierr = NF_INQ_VARID(ncid,'RUGOS',varid)
    259       if (ierr.ne.0) then
    260         print *, NF_STRERROR(ierr)
    261         STOP
    262       ENDIF
    263       ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
    264       if (ierr.ne.0) then
    265         print *, NF_STRERROR(ierr)
    266         STOP
    267       ENDIF
    268       ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
    269       if (ierr.ne.0) then
    270         print *, NF_STRERROR(ierr)
    271         STOP
    272       ENDIF
    273       print*,'variable ', namedim, 'dimension ', imdep
    274       ierr = NF_INQ_VARID(ncid,namedim,dimid)
    275       if (ierr.ne.0) then
    276         print *, NF_STRERROR(ierr)
    277         STOP
    278       ENDIF
    279 
    280       ALLOCATE( dlon_ini(imdep) )
    281       ALLOCATE(     dlon(imdep) )
    282 
    283 #ifdef NC_DOUBLE
    284       ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini)
    285 #else
    286       ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini)
    287 #endif
    288       if (ierr.ne.0) then
    289         print *, NF_STRERROR(ierr)
    290         STOP
    291       ENDIF
    292       ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
    293       if (ierr.ne.0) then
    294         print *, NF_STRERROR(ierr)
    295         STOP
    296       ENDIF
    297       print*,'variable ', namedim, 'dimension ', jmdep
    298       ierr = NF_INQ_VARID(ncid,namedim,dimid)
    299       if (ierr.ne.0) then
    300         print *, NF_STRERROR(ierr)
    301         STOP
    302       ENDIF
    303 
    304       ALLOCATE( dlat_ini(jmdep) )
    305       ALLOCATE(     dlat(jmdep) )
    306 
    307 #ifdef NC_DOUBLE
    308       ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini)
    309 #else
    310       ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini)
    311 #endif
    312       if (ierr.ne.0) then
    313         print *, NF_STRERROR(ierr)
    314         STOP
    315       ENDIF
    316       ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
    317       if (ierr.ne.0) then
    318         print *, NF_STRERROR(ierr)
    319         STOP
    320       ENDIF
    321       print*,'variable ', namedim, 'dimension ', lmdep
    322       ierr = NF_INQ_VARID(ncid,namedim,dimid)
    323       if (ierr.ne.0) then
    324         print *, NF_STRERROR(ierr)
    325         STOP
    326       ENDIF
    327 #ifdef NC_DOUBLE
    328       ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
    329 #else
    330       ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
    331 #endif
    332       if (ierr.ne.0) then
    333         print *, NF_STRERROR(ierr)
    334         STOP
    335       ENDIF
    336 c
    337       ALLOCATE( champ(imdep*jmdep) )
    338 
    339       DO  200 l = 1, lmdep
    340          dimfirst(1) = 1
    341          dimfirst(2) = 1
    342          dimfirst(3) = l
    343 c
    344          dimlast(1) = imdep
    345          dimlast(2) = jmdep
    346          dimlast(3) = 1
    347 c
    348          PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
    349          print*,dimfirst,dimlast
    350 #ifdef NC_DOUBLE
    351          ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
    352 #else
    353          ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
    354 #endif
    355          if (ierr.ne.0) then
    356            print *, NF_STRERROR(ierr)
    357            STOP
    358          ENDIF
    359    
    360         title = 'Rugosite Amip '
    361 c
    362         CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
    363      .                      dlon, dlat, champ, interbar          )
    364 
    365        IF ( interbar )   THEN
    366          DO j = 1, imdep * jmdep
    367            champ(j) = LOG(champ(j))
    368          ENDDO
    369 
    370         IF( l.EQ.1 )  THEN
    371          WRITE(6,*) '-------------------------------------------------',
    372      ,'------------------------'
    373          WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
    374      , ' pour la rugosite $$$ '
    375          WRITE(6,*) '-------------------------------------------------',
    376      ,'------------------------'
    377         ENDIF
    378         CALL inter_barxy ( imdep,jmdep -1,dlon,dlat,champ ,
    379      ,                  iim,jjm,rlonu,rlatv, jjp1,champint )
    380          DO j=1,jjp1
    381           DO i=1,iim
    382            champint(i,j)=EXP(champint(i,j))
    383           ENDDO
    384          ENDDO
    385 
    386          DO j = 1, jjp1
    387            DO i = 1, iim
    388              IF(NINT(mask(i,j)).NE.1)  THEN
    389                champint( i,j ) = 0.001
    390              ENDIF
    391            ENDDO
    392          ENDDO
    393       ELSE
    394          CALL rugosite(imdep, jmdep, dlon, dlat, champ,
    395      .             iim, jjp1, rlonv, rlatu, champint, mask)
    396       ENDIF
    397          DO j = 1,jjp1
    398          DO i = 1, iim
    399             champtime (i,j,l) = champint(i,j)
    400          ENDDO
    401          ENDDO
    402 200      CONTINUE
    403 c
    404       DO l = 1, lmdep
    405          timeyear(l) = timecoord(l)
    406       ENDDO
    407 
    408       PRINT 222, timeyear(:lmdep)
    409 222   FORMAT(2x,' Time year ',10f6.1)
    410 c
    411        
    412       PRINT*, 'Interpolation temporelle dans l annee'
    413 
    414       DO j = 1, jjp1
    415       DO i = 1, iim
    416           DO l = 1, lmdep
    417             ax(l) = timeyear(l)
    418             ay(l) = champtime (i,j,l)
    419           ENDDO
    420           CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
    421           DO k = 1, 360
    422             time = FLOAT(k-1)
    423             CALL SPLINT(ax,ay,yder,lmdep,time,by)
    424             champan(i,j,k) = by
    425           ENDDO
    426       ENDDO
    427       ENDDO
    428       DO k = 1, 360
    429       DO j = 1, jjp1
    430          champan(iip1,j,k) = champan(1,j,k)
    431       ENDDO
    432         IF ( k.EQ.10 )  THEN
    433           DO j = 1, jjp1
    434             CALL minmax( iip1,champan(1,j,10),chmin,chmax )
    435             PRINT *,' Rugosite au temps 10 ', chmin,chmax,j
    436           ENDDO
    437         ENDIF
    438       ENDDO
    439 c
    440       DO k = 1, 360
    441          CALL gr_dyn_fi(1,iip1,jjp1,klon,champan(1,1,k), phy_rug(1,k))
    442       ENDDO
    443 c
    444       ierr = NF_CLOSE(ncid)
    445 
    446        DEALLOCATE( dlon      )
    447        DEALLOCATE( dlon_ini  )
    448        DEALLOCATE( dlat      )
    449        DEALLOCATE( dlat_ini  )
    450        DEALLOCATE( champ     )
    451 c
    452 c
    453 C Traitement de la glace oceanique
    454 c
    455       PRINT*, 'Traitement de la glace oceanique'
    456 
    457       ierr = NF_OPEN('amipbc_sic_1x1.nc', NF_NOWRITE, ncid)
    458       if (ierr.ne.0) THEN
    459         ierr = NF_OPEN('amipbc_sic_1x1_clim.nc', NF_NOWRITE, ncid)
    460         if (ierr.ne.0) THEN
    461           print *, NF_STRERROR(ierr)
    462           STOP
    463         endif
    464       ENDIF
    465 
    466 cIM22/02/2002
    467 cIM07/03/2002 AMIP.nc & amip79to95.nc
    468 cIM   ierr = NF_INQ_VARID(ncid,'SEA_ICE',varid)
    469 cIM07/03/2002 amipbc_sic_1x1_clim.nc & amipbc_sic_1x1.nc
    470       ierr = NF_INQ_VARID(ncid,'sicbcs',varid)
    471 cIM22/02/2002
    472       if (ierr.ne.0) then
    473         print *, NF_STRERROR(ierr),'sicbcs'
    474         STOP
    475       ENDIF
    476       ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
    477       if (ierr.ne.0) then
    478         print *, NF_STRERROR(ierr)
    479         STOP
    480       ENDIF
    481       ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
    482       if (ierr.ne.0) then
    483         print *, NF_STRERROR(ierr)
    484         STOP
    485       ENDIF
    486       print*,'variable ', namedim, 'dimension ', imdep
    487       ierr = NF_INQ_VARID(ncid,namedim,dimid)
    488       if (ierr.ne.0) then
    489         print *, NF_STRERROR(ierr)
    490         STOP
    491       ENDIF
    492 
    493       ALLOCATE ( dlon_ini(imdep) )
    494       ALLOCATE (     dlon(imdep) )
    495 
    496 #ifdef NC_DOUBLE
    497       ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini)
    498 #else
    499       ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini)
    500 #endif
    501       if (ierr.ne.0) then
    502         print *, NF_STRERROR(ierr)
    503         STOP
    504       ENDIF
    505       ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
    506       if (ierr.ne.0) then
    507         print *, NF_STRERROR(ierr)
    508         STOP
    509       ENDIF
    510       print*,'variable ', namedim, jmdep
    511       ierr = NF_INQ_VARID(ncid,namedim,dimid)
    512       if (ierr.ne.0) then
    513         print *, NF_STRERROR(ierr)
    514         STOP
    515       ENDIF
    516 
    517       ALLOCATE ( dlat_ini(jmdep) )
    518       ALLOCATE (     dlat(jmdep) )
    519 
    520 #ifdef NC_DOUBLE
    521       ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini)
    522 #else
    523       ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini)
    524 #endif
    525       if (ierr.ne.0) then
    526         print *, NF_STRERROR(ierr)
    527         STOP
    528       ENDIF
    529       ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
    530       if (ierr.ne.0) then
    531         print *, NF_STRERROR(ierr)
    532         STOP
    533       ENDIF
    534       print*,'variable ', namedim, lmdep
    535 cIM28/02/2002
    536 cPM28/02/2002 : nouvelle coord temporelle fichiers AMIP pas en jours
    537 c               Ici on suppose qu'on a 12 mois (de 30 jours).
    538       IF (lmdep.NE.12) THEN
    539           print *, 'Unknown AMIP file: not 12 months ?'
    540           STOP
    541        ENDIF
    542 cIM28/02/2002
    543       ierr = NF_INQ_VARID(ncid,namedim,dimid)
    544       if (ierr.ne.0) then
    545         print *, NF_STRERROR(ierr)
    546         STOP
    547       ENDIF
    548 #ifdef NC_DOUBLE
    549       ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
    550 #else
    551       ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
    552 #endif
    553       if (ierr.ne.0) then
    554         print *, NF_STRERROR(ierr)
    555         STOP
    556       ENDIF
    557 c
    558       ALLOCATE ( champ(imdep*jmdep) )
    559 
    560       DO l = 1, lmdep
    561          dimfirst(1) = 1
    562          dimfirst(2) = 1
    563          dimfirst(3) = l
    564 c
    565          dimlast(1) = imdep
    566          dimlast(2) = jmdep
    567          dimlast(3) = 1
    568 c
    569          PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
    570 #ifdef NC_DOUBLE
    571          ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
    572 #else
    573          ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
    574 #endif
    575          if (ierr.ne.0) then
    576            print *, NF_STRERROR(ierr)
    577            STOP
    578          ENDIF
    579  
    580          title = 'Sea-ice Amip '
    581 c
    582          CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
    583      .                        dlon, dlat, champ, interbar          )
    584 c
    585       IF( oldice )  THEN
    586                  CALL sea_ice(imdep, jmdep, dlon, dlat, champ,
    587      .             iim, jjp1, rlonv, rlatu, champint )
    588       ELSEIF ( interbar )  THEN
    589        IF( l.EQ.1 )  THEN
    590         WRITE(6,*) '-------------------------------------------------',
    591      ,'------------------------'
    592         WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
    593      , ' pour Sea-ice Amip  $$$ '
    594         WRITE(6,*) '-------------------------------------------------',
    595      ,'------------------------'
    596        ENDIF
    597 cIM07/03/2002
    598 cIM22/02/2002 : Sea-ice Amip entre 0. et 1.
    599 cIM    PRINT*,'SUB. limit_netcdf.F IM : Sea-ice et SST Amip_new clim'
    600 cIM   DO j = 1, imdep * jmdep
    601 cIM28/02/2002 <==PM         champ(j) = champ(j)/100.
    602 cIM14/03/2002      champ(j) = max(0.0,(min(1.0, (champ(j)/100.) )))
    603 cIM      champ(j) = amax1(0.0,(amin1(1.0, (champ(j)/100.) )))
    604 cIM   ENDDO
    605 cIM22/02/2002
    606          CALL inter_barxy ( imdep,jmdep -1,dlon, dlat ,
    607      ,     champ, iim, jjm, rlonu, rlatv, jjp1, champint )
    608       ELSE
    609          CALL sea_ice(imdep, jmdep, dlon, dlat, champ,
    610      .             iim, jjp1, rlonv, rlatu, champint )
    611       ENDIF
    612          DO j = 1,jjp1
    613          DO i = 1, iim
    614             champtime (i,j,l) = champint(i,j)
    615          ENDDO
    616          ENDDO
    617       ENDDO
    618 c
    619       DO l = 1, lmdep
    620 cIM28/02/2002 <== PM  timeyear(l) = timecoord(l)
    621 cIM      timeyear(l) = timecoord(l)
    622 cIM07/03/2002     
    623          timeyear(l) = tmidmonth(l)
    624       ENDDO
    625       PRINT 222,  timeyear(:lmdep)
    626 c
    627       PRINT*, 'Interpolation temporelle'
    628       DO j = 1, jjp1
    629       DO i = 1, iim
    630           DO l = 1, lmdep
    631             ax(l) = timeyear(l)
    632             ay(l) = champtime (i,j,l)
    633           ENDDO
    634           CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
    635           DO k = 1, 360
    636             time = FLOAT(k-1)
    637             CALL SPLINT(ax,ay,yder,lmdep,time,by)
    638             champan(i,j,k) = by
    639           ENDDO
    640       ENDDO
    641       ENDDO
    642       DO k = 1, 360
    643       DO j = 1, jjp1
    644          champan(iip1, j, k) = champan(1, j, k)
    645       ENDDO
    646         IF ( k.EQ.10 )  THEN
    647           DO j = 1, jjp1
    648             CALL minmax( iip1,champan(1,j,10),chmin,chmax )
    649             PRINT *,' Sea ice au temps 10 ', chmin,chmax,j
    650           ENDDO
    651         ENDIF
    652       ENDDO
    653 c
    654 cIM14/03/2002 : Sea-ice Amip entre 0. et 1.
    655       PRINT*,'SUB. limit_netcdf.F IM : Sea-ice Amipbc '
    656       DO k = 1, 360
    657       DO j = 1, jjp1
    658       DO i = 1, iim
    659         champan(i, j, k) =
    660      $ amax1(0.0,(amin1(1.0,(champan(i, j, k)/100.))))
    661       ENDDO
    662         champan(iip1, j, k) = champan(1, j, k)
    663       ENDDO
    664       ENDDO
    665 cIM14/03/2002
    666 
    667       DO k = 1, 360
    668          CALL gr_dyn_fi(1, iip1, jjp1, klon,
    669      .                  champan(1,1,k), phy_ice)
    670         IF ( newlmt) THEN
    671 
    672 CPB  en attendant de mettre fraction de terre
    673 c
    674           WHERE(phy_ice(1:klon) .GE. 1.) phy_ice(1 : klon) = 1.
    675           WHERE(phy_ice(1:klon) .LT. EPSFRA) phy_ice(1 : klon) = 0.
    676 c
    677           IF (fracterre ) THEN
    678 c            WRITE(*,*) 'passe dans cas fracterre'
    679             pctsrf_t(:,is_ter,k) = pctsrf(:,is_ter)
    680             pctsrf_t(:,is_lic,k) = pctsrf(:,is_lic)
    681             pctsrf_t(1:klon,is_sic,k) =   phy_ice(1:klon)
    682      $            - pctsrf_t(1:klon,is_lic,k)
    683 c Il y a des cas ou il y a de la glace dans landiceref et pas dans AMIP
    684             WHERE (pctsrf_t(1:klon,is_sic,k) .LE. 0)
    685               pctsrf_t(1:klon,is_sic,k) = 0.
    686             END WHERE
    687             WHERE( 1. - zmasq(1:klon) .LT. EPSFRA)
    688               pctsrf_t(1:klon,is_sic,k) = 0.
    689               pctsrf_t(1:klon,is_oce,k) = 0.
    690             END WHERE
    691             DO i = 1, klon
    692               IF ( 1. - zmasq(i) .GT. EPSFRA) THEN
    693                 IF ( pctsrf_t(i,is_sic,k) .GE. 1 - zmasq(i)) THEN
    694                   pctsrf_t(i,is_sic,k) = 1 - zmasq(i)
    695                   pctsrf_t(i,is_oce,k) = 0.
    696                 ELSE
    697                   pctsrf_t(i,is_oce,k) = 1 - zmasq(i)
    698      $                    - pctsrf_t(i,is_sic,k)
    699                   IF (pctsrf_t(i,is_oce,k) .LT. EPSFRA) THEN
    700                     pctsrf_t(i,is_oce,k) = 0.
    701                     pctsrf_t(i,is_sic,k) = 1 - zmasq(i)
    702                   ENDIF
    703                 ENDIF
    704               ENDIF 
    705               if (pctsrf_t(i,is_oce,k) .lt. 0.) then
    706                 WRITE(*,*) 'pb sous maille au point : i,k '
    707      $              , i,k,pctsrf_t(:,is_oce,k)
    708               ENDIF
    709               IF ( abs( pctsrf_t(i, is_ter,k) + pctsrf_t(i, is_lic,k) +
    710      $          pctsrf_t(i, is_oce,k) + pctsrf_t(i, is_sic,k)  - 1.)
    711      $            .GT. EPSFRA) THEN
    712                   WRITE(*,*) 'physiq : pb sous surface au point ', i,
    713      $                pctsrf_t(i, 1 : nbsrf,k), phy_ice(i)
    714               ENDIF
    715             END DO
    716           ELSE
    717             DO i = 1, klon
    718               pctsrf_t(i,is_ter,k) = pctsrf(i,is_ter)
    719               IF (NINT(pctsrf(i,is_ter)).EQ.1 ) THEN
    720                 pctsrf_t(i,is_sic,k) = 0.
    721                 pctsrf_t(i,is_oce,k) = 0.                 
    722                 IF(phy_ice(i) .GE. 1.e-5) THEN
    723                   pctsrf_t(i,is_lic,k) = phy_ice(i)
    724                   pctsrf_t(i,is_ter,k) = pctsrf_t(i,is_ter,k)
    725      .                                   - pctsrf_t(i,is_lic,k)
    726                 ELSE
    727                   pctsrf_t(i,is_lic,k) = 0.
    728                 ENDIF
    729               ELSE
    730                 pctsrf_t(i,is_lic,k) = 0.
    731                 IF(phy_ice(i) .GE. 1.e-5) THEN
    732                   pctsrf_t(i,is_ter,k) = 0.
    733                   pctsrf_t(i,is_sic,k) = phy_ice(i)
    734                   pctsrf_t(i,is_oce,k) = 1. - pctsrf_t(i,is_sic,k)
    735                 ELSE
    736                   pctsrf_t(i,is_sic,k) = 0.
    737                   pctsrf_t(i,is_oce,k) = 1.                     
    738                 ENDIF
    739               ENDIF
    740               verif = pctsrf_t(i,is_ter,k) +
    741      .                pctsrf_t(i,is_oce,k) +
    742      .                pctsrf_t(i,is_sic,k) +
    743      .                pctsrf_t(i,is_lic,k)
    744               IF ( verif .LT. 1. - 1.e-5 .OR.
    745      $             verif .GT. 1 + 1.e-5) THEN 
    746                 WRITE(*,*) 'pb sous maille au point : i,k,verif '
    747      $                    , i,k,verif
    748               ENDIF
    749             END DO
    750           ENDIF
    751         ELSE 
    752           DO i = 1, klon
    753             phy_nat(i,k) = phy_nat0(i)
    754             IF ( (phy_ice(i) - 0.5).GE.1.e-5 ) THEN
    755               IF (NINT(phy_nat0(i)).EQ.0) THEN
    756                 phy_nat(i,k) = 3.0
    757               ELSE
    758                 phy_nat(i,k) = 2.0
    759               ENDIF
    760             ENDIF
    761             IF( NINT(phy_nat(i,k)).EQ.0 ) THEN
    762               IF ( phy_rug(i,k).NE.0.001 ) phy_rug(i,k) = 0.001
    763             ENDIF
    764           END DO
    765         ENDIF
    766       ENDDO
    767 c
    768 
    769       ierr = NF_CLOSE(ncid)
    770 c
    771        DEALLOCATE( dlon      )
    772        DEALLOCATE( dlon_ini  )
    773        DEALLOCATE( dlat      )
    774        DEALLOCATE( dlat_ini  )
    775        DEALLOCATE( champ     )
    776 
    777 477    continue
    778 c
    779 C Traitement de la sst
    780 c
    781       PRINT*, 'Traitement de la sst'
    782 c     ierr = NF_OPEN('AMIP_SST.nc', NF_NOWRITE, ncid)
    783       ierr = NF_OPEN('amipbc_sst_1x1.nc', NF_NOWRITE, ncid)
    784       if (ierr.ne.0) THEN
    785         ierr = NF_OPEN('amipbc_sst_1x1_clim.nc', NF_NOWRITE, ncid)
    786         if (ierr.ne.0) THEN
    787           print *, NF_STRERROR(ierr)
    788           STOP
    789         endif
    790       ENDIF
    791 
    792 cIM22/02/2002
    793 cIM   ierr = NF_INQ_VARID(ncid,'SST',varid)
    794       ierr = NF_INQ_VARID(ncid,'tosbcs',varid)
    795 cIM22/02/2002
    796       if (ierr.ne.0) then
    797         print *, NF_STRERROR(ierr)
    798         STOP
    799       ENDIF
    800       ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
    801       if (ierr.ne.0) then
    802         print *, NF_STRERROR(ierr)
    803         STOP
    804       ENDIF
    805       ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
    806       if (ierr.ne.0) then
    807         print *, NF_STRERROR(ierr)
    808         STOP
    809       ENDIF
    810       print*,'variable SST ', namedim,'dimension ', imdep
    811       ierr = NF_INQ_VARID(ncid,namedim,dimid)
    812       if (ierr.ne.0) then
    813         print *, NF_STRERROR(ierr)
    814         STOP
    815       ENDIF
    816  
    817       ALLOCATE( dlon_ini(imdep) )
    818       ALLOCATE(     dlon(imdep) )
    819 
    820 #ifdef NC_DOUBLE
    821       ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini)
    822 #else
    823       ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini)
    824 #endif
    825 
    826       if (ierr.ne.0) then
    827         print *, NF_STRERROR(ierr)
    828         STOP
    829       ENDIF
    830       ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
    831       if (ierr.ne.0) then
    832         print *, NF_STRERROR(ierr)
    833         STOP
    834       ENDIF
    835       print*,'variable SST ', namedim, 'dimension ', jmdep
    836       ierr = NF_INQ_VARID(ncid,namedim,dimid)
    837       if (ierr.ne.0) then
    838         print *, NF_STRERROR(ierr)
    839         STOP
    840       ENDIF
    841 
    842       ALLOCATE( dlat_ini(jmdep) )
    843       ALLOCATE(     dlat(jmdep) )
    844 
    845 #ifdef NC_DOUBLE
    846       ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini)
    847 #else
    848       ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini)
    849 #endif
    850       if (ierr.ne.0) then
    851         print *, NF_STRERROR(ierr)
    852         STOP
    853       ENDIF
    854       ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
    855       if (ierr.ne.0) then
    856         print *, NF_STRERROR(ierr)
    857         STOP
    858       ENDIF
    859       print*,'variable ', namedim, 'dimension ', lmdep
    860 cIM28/02/2002
    861 cPM28/02/2002 : nouvelle coord temporelle fichiers AMIP pas en jours
    862 c               Ici on suppose qu'on a 12 mois (de 30 jours).
    863       IF (lmdep.NE.12) THEN
    864           print *, 'Unknown AMIP file: not 12 months ?'
    865           STOP
    866        ENDIF
    867 cIM28/02/2002
    868       ierr = NF_INQ_VARID(ncid,namedim,dimid)
    869       if (ierr.ne.0) then
    870         print *, NF_STRERROR(ierr)
    871         STOP
    872       ENDIF
    873 #ifdef NC_DOUBLE
    874       ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
    875 #else
    876       ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
    877 #endif
    878       if (ierr.ne.0) then
    879         print *, NF_STRERROR(ierr)
    880         STOP
    881       ENDIF
    882 
    883        ALLOCATE( champ(imdep*jmdep) )
    884        IF( extrap )   THEN
    885          ALLOCATE ( work(imdep,jmdep) )
    886        ENDIF
    887 c
    888       DO l = 1, lmdep
    889          dimfirst(1) = 1
    890          dimfirst(2) = 1
    891          dimfirst(3) = l
    892 c
    893          dimlast(1) = imdep
    894          dimlast(2) = jmdep
    895          dimlast(3) = 1
    896 c
    897          PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
    898 #ifdef NC_DOUBLE
    899          ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
    900 #else
    901          ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
    902 #endif
    903          if (ierr.ne.0) then
    904            print *, NF_STRERROR(ierr)
    905            STOP
    906          ENDIF
    907 
    908          title='Sst Amip'
    909 c
    910          CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
    911      .                            dlon, dlat, champ, interbar     )
    912        IF ( extrap )  THEN
    913         CALL extrapol(champ, imdep, jmdep, 999999.,.TRUE.,.TRUE.,2,work)
    914        ENDIF
    915 c
    916 
    917       IF ( interbar )  THEN
    918         IF( l.EQ.1 )  THEN
    919          WRITE(6,*) '-------------------------------------------------',
    920      ,'------------------------'
    921          WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
    922      , ' pour la Sst Amip $$$ '
    923          WRITE(6,*) '-------------------------------------------------',
    924      ,'------------------------'
    925         ENDIF
    926        CALL inter_barxy ( imdep,jmdep -1,dlon, dlat ,
    927      , champ, iim, jjm, rlonu, rlatv, jjp1, champint )
    928       ELSE
    929        CALL grille_m(imdep, jmdep, dlon, dlat, champ,
    930      .          iim, jjp1, rlonv, rlatu, champint   )
    931       ENDIF
    932 
    933          DO j = 1,jjp1
    934          DO i = 1, iim
    935             champtime (i,j,l) = champint(i,j)
    936          ENDDO
    937          ENDDO
    938       ENDDO
    939 c
    940       DO l = 1, lmdep
    941 cIM28/02/2002 <==PM  timeyear(l) = timecoord(l)
    942          timeyear(l) = tmidmonth(l)
    943       ENDDO
    944       print 222,  timeyear(:lmdep)
    945 c
    946 C interpolation temporelle
    947       DO j = 1, jjp1
    948       DO i = 1, iim
    949           DO l = 1, lmdep
    950             ax(l) = timeyear(l)
    951             ay(l) = champtime (i,j,l)
    952           ENDDO
    953           CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
    954           DO k = 1, 360
    955             time = FLOAT(k-1)
    956             CALL SPLINT(ax,ay,yder,lmdep,time,by)
    957             champan(i,j,k) = by
    958           ENDDO
    959       ENDDO
    960       ENDDO
    961       DO k = 1, 360
    962       DO j = 1, jjp1
    963          champan(iip1,j,k) = champan(1,j,k)
    964       ENDDO
    965         IF ( k.EQ.10 )  THEN
    966           DO j = 1, jjp1
    967             CALL minmax( iip1,champan(1,j,10),chmin,chmax )
    968             PRINT *,' SST au temps 10 ', chmin,chmax,j
    969           ENDDO
    970         ENDIF
    971       ENDDO
    972 c
    973 cIM14/03/2002 : SST amipbc greater then 271.38
    974       PRINT*,'SUB. limit_netcdf.F IM : SST Amipbc >= 271.38 '
    975       DO k = 1, 360
    976       DO j = 1, jjp1
    977       DO i = 1, iim
    978          champan(i, j, k) = amax1(champan(i, j, k), 271.38)
    979       ENDDO
    980          champan(iip1, j, k) = champan(1, j, k)
    981       ENDDO
    982       ENDDO
    983 cIM14/03/2002
    984       DO k = 1, 360
    985          CALL gr_dyn_fi(1, iip1, jjp1, klon,
    986      .                  champan(1,1,k), phy_sst(1,k))
    987       ENDDO
    988 c
    989       ierr = NF_CLOSE(ncid)
    990 c
    991        DEALLOCATE( dlon      )
    992        DEALLOCATE( dlon_ini  )
    993        DEALLOCATE( dlat      )
    994        DEALLOCATE( dlat_ini  )
    995        DEALLOCATE( champ     )
    996 c
    997 C Traitement de l'albedo
    998 c
    999       PRINT*, 'Traitement de l albedo'
    1000       ierr = NF_OPEN('Albedo.nc', NF_NOWRITE, ncid)
    1001       if (ierr.ne.0) then
    1002         print *, NF_STRERROR(ierr)
    1003         STOP
    1004       ENDIF
    1005       ierr = NF_INQ_VARID(ncid,'ALBEDO',varid)
    1006       if (ierr.ne.0) then
    1007         print *, NF_STRERROR(ierr)
    1008         STOP
    1009       ENDIF
    1010       ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
    1011       if (ierr.ne.0) then
    1012         print *, NF_STRERROR(ierr)
    1013         STOP
    1014       ENDIF
    1015       ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
    1016       if (ierr.ne.0) then
    1017         print *, NF_STRERROR(ierr)
    1018         STOP
    1019       ENDIF
    1020       print*,'variable ', namedim, 'dimension ', imdep
    1021       ierr = NF_INQ_VARID(ncid,namedim,dimid)
    1022       if (ierr.ne.0) then
    1023         print *, NF_STRERROR(ierr)
    1024         STOP
    1025       ENDIF
    1026 
    1027       ALLOCATE ( dlon_ini(imdep) )
    1028       ALLOCATE (     dlon(imdep) )
    1029 
    1030 #ifdef NC_DOUBLE
    1031       ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini)
    1032 #else
    1033       ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini)
    1034 #endif
    1035       if (ierr.ne.0) then
    1036         print *, NF_STRERROR(ierr)
    1037         STOP
    1038       ENDIF
    1039       ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
    1040       if (ierr.ne.0) then
    1041         print *, NF_STRERROR(ierr)
    1042         STOP
    1043       ENDIF
    1044       print*,'variable ', namedim, 'dimension ', jmdep
    1045       ierr = NF_INQ_VARID(ncid,namedim,dimid)
    1046       if (ierr.ne.0) then
    1047         print *, NF_STRERROR(ierr)
    1048         STOP
    1049       ENDIF
    1050 
    1051       ALLOCATE ( dlat_ini(jmdep) )
    1052       ALLOCATE (     dlat(jmdep) )
    1053 
    1054 #ifdef NC_DOUBLE
    1055       ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini)
    1056 #else
    1057       ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini)
    1058 #endif
    1059       if (ierr.ne.0) then
    1060         print *, NF_STRERROR(ierr)
    1061         STOP
    1062       ENDIF
    1063       ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
    1064       if (ierr.ne.0) then
    1065         print *, NF_STRERROR(ierr)
    1066         STOP
    1067       ENDIF
    1068       print*,'variable ', namedim, 'dimension ', lmdep
    1069       ierr = NF_INQ_VARID(ncid,namedim,dimid)
    1070       if (ierr.ne.0) then
    1071         print *, NF_STRERROR(ierr)
    1072         STOP
    1073       ENDIF
    1074 #ifdef NC_DOUBLE
    1075       ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
    1076 #else
    1077       ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
    1078 #endif
    1079       if (ierr.ne.0) then
    1080         print *, NF_STRERROR(ierr)
    1081         STOP
    1082       ENDIF
    1083 c
    1084       ALLOCATE ( champ(imdep*jmdep) )
    1085 
    1086       DO l = 1, lmdep
    1087          dimfirst(1) = 1
    1088          dimfirst(2) = 1
    1089          dimfirst(3) = l
    1090 c
    1091          dimlast(1) = imdep
    1092          dimlast(2) = jmdep
    1093          dimlast(3) = 1
    1094 c
    1095          PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
    1096 #ifdef NC_DOUBLE
    1097          ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
    1098 #else
    1099          ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
    1100 #endif
    1101          if (ierr.ne.0) then
    1102            print *, NF_STRERROR(ierr)
    1103            STOP
    1104          ENDIF
    1105 
    1106          title='Albedo Amip'
    1107 c
    1108          CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
    1109      .                            dlon, dlat, champ, interbar      )
    1110 c
    1111 c
    1112       IF ( interbar )  THEN
    1113         IF( l.EQ.1 )  THEN
    1114          WRITE(6,*) '-------------------------------------------------',
    1115      ,'------------------------'
    1116          WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
    1117      , ' pour l Albedo Amip $$$ '
    1118          WRITE(6,*) '-------------------------------------------------',
    1119      ,'------------------------'
    1120         ENDIF
    1121 
    1122        CALL inter_barxy ( imdep,jmdep -1,dlon, dlat ,
    1123      , champ, iim, jjm, rlonu, rlatv, jjp1, champint )
    1124       ELSE
    1125        CALL grille_m(imdep, jmdep, dlon, dlat, champ,
    1126      .          iim, jjp1, rlonv, rlatu, champint   )
    1127       ENDIF
    1128 c
    1129          DO j = 1,jjp1
    1130          DO i = 1, iim
    1131             champtime (i, j, l) = champint(i, j)
    1132          ENDDO
    1133          ENDDO
    1134       ENDDO
    1135 c
    1136       DO l = 1, lmdep
    1137          timeyear(l) = timecoord(l)
    1138       ENDDO
    1139       print 222,  timeyear(:lmdep)
    1140 c
    1141 C interpolation temporelle
    1142       DO j = 1, jjp1
    1143       DO i = 1, iim
    1144           DO l = 1, lmdep
    1145             ax(l) = timeyear(l)
    1146             ay(l) = champtime (i, j, l)
    1147           ENDDO
    1148           CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
    1149           DO k = 1, 360
    1150             time = FLOAT(k-1)
    1151             CALL SPLINT(ax,ay,yder,lmdep,time,by)
    1152             champan(i,j,k) = by
    1153           ENDDO
    1154       ENDDO
    1155       ENDDO
    1156       DO k = 1, 360
    1157       DO j = 1, jjp1
    1158          champan(iip1, j, k) = champan(1, j, k)
    1159       ENDDO
    1160         IF ( k.EQ.10 )  THEN
    1161           DO j = 1, jjp1
    1162             CALL minmax( iip1,champan(1,j,10),chmin,chmax )
    1163             PRINT *,' Albedo au temps 10 ', chmin,chmax,j
    1164           ENDDO
    1165         ENDIF
    1166       ENDDO
    1167 c
    1168       DO k = 1, 360
    1169          CALL gr_dyn_fi(1, iip1, jjp1, klon,
    1170      .                  champan(1,1,k), phy_alb(1,k))
    1171       ENDDO
    1172 c
    1173       ierr = NF_CLOSE(ncid)
    1174 c
    1175 c
    1176       DO k = 1, 360
    1177       DO i = 1, klon
    1178          phy_bil(i,k) = 0.0
    1179       ENDDO
    1180       ENDDO
    1181 c
    1182       PRINT*, 'Ecriture du fichier limit'
    1183 c
    1184       ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid)
    1185 c
    1186       ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30,
    1187      .                       "Fichier conditions aux limites")
    1188       ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
    1189       ierr = NF_DEF_DIM (nid, "time", NF_UNLIMITED, ntim)
    1190 c
    1191       dims(1) = ndim
    1192       dims(2) = ntim
    1193 c
    1194 #ifdef NC_DOUBLE
    1195       ierr = NF_DEF_VAR (nid, "TEMPS", NF_DOUBLE, 1,ntim, id_tim)
    1196 #else
    1197       ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim)
    1198 #endif
    1199       ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17,
    1200      .                        "Jour dans l annee")
    1201       IF (newlmt) THEN
    1202 c
    1203 #ifdef NC_DOUBLE
    1204         ierr = NF_DEF_VAR (nid, "FOCE", NF_DOUBLE, 2,dims, id_FOCE)
    1205 #else
    1206         ierr = NF_DEF_VAR (nid, "FOCE", NF_FLOAT, 2,dims, id_FOCE)
    1207 #endif
    1208         ierr = NF_PUT_ATT_TEXT (nid, id_FOCE, "title", 14,
    1209      .                      "Fraction ocean")
    1210 c
    1211 #ifdef NC_DOUBLE
    1212         ierr = NF_DEF_VAR (nid, "FSIC", NF_DOUBLE, 2,dims, id_FSIC)
    1213 #else
    1214         ierr = NF_DEF_VAR (nid, "FSIC", NF_FLOAT, 2,dims, id_FSIC)
    1215 #endif
    1216         ierr = NF_PUT_ATT_TEXT (nid, id_FSIC, "title", 21,
    1217      .                      "Fraction glace de mer")
    1218 c
    1219 #ifdef NC_DOUBLE
    1220         ierr = NF_DEF_VAR (nid, "FTER", NF_DOUBLE, 2,dims, id_FTER)
    1221 #else
    1222         ierr = NF_DEF_VAR (nid, "FTER", NF_FLOAT, 2,dims, id_FTER)
    1223 #endif
    1224         ierr = NF_PUT_ATT_TEXT (nid, id_FTER, "title", 14,
    1225      .                      "Fraction terre")
    1226 c
    1227 #ifdef NC_DOUBLE
    1228         ierr = NF_DEF_VAR (nid, "FLIC", NF_DOUBLE, 2,dims, id_FLIC)
    1229 #else
    1230         ierr = NF_DEF_VAR (nid, "FLIC", NF_FLOAT, 2,dims, id_FLIC)
    1231 #endif
    1232         ierr = NF_PUT_ATT_TEXT (nid, id_FLIC, "title", 17,
    1233      .                      "Fraction land ice")
    1234 c
    1235       ELSE
    1236 #ifdef NC_DOUBLE
    1237         ierr = NF_DEF_VAR (nid, "NAT", NF_DOUBLE, 2,dims, id_NAT)
    1238 #else
    1239         ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT)
    1240 #endif
    1241         ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23,
    1242      .                      "Nature du sol (0,1,2,3)")
    1243       ENDIF
    1244 #ifdef NC_DOUBLE
    1245       ierr = NF_DEF_VAR (nid, "SST", NF_DOUBLE, 2,dims, id_SST)
    1246 #else
    1247       ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST)
    1248 #endif
    1249       ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35,
    1250      .                      "Temperature superficielle de la mer")
    1251 #ifdef NC_DOUBLE
    1252       ierr = NF_DEF_VAR (nid, "BILS", NF_DOUBLE, 2,dims, id_BILS)
    1253 #else
    1254       ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS)
    1255 #endif
    1256       ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32,
    1257      .                        "Reference flux de chaleur au sol")
    1258 #ifdef NC_DOUBLE
    1259       ierr = NF_DEF_VAR (nid, "ALB", NF_DOUBLE, 2,dims, id_ALB)
    1260 #else
    1261       ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB)
    1262 #endif
    1263       ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19,
    1264      .                        "Albedo a la surface")
    1265 #ifdef NC_DOUBLE
    1266       ierr = NF_DEF_VAR (nid, "RUG", NF_DOUBLE, 2,dims, id_RUG)
    1267 #else
    1268       ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG)
    1269 #endif
    1270       ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8,
    1271      .                        "Rugosite")
    1272 c
    1273       ierr = NF_ENDDEF(nid)
    1274 c
    1275       DO k = 1, 360
    1276 c
    1277       debut(1) = 1
    1278       debut(2) = k
    1279       epais(1) = klon
    1280       epais(2) = 1
    1281 c
    1282 #ifdef NC_DOUBLE
    1283       ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k))
    1284 c
    1285       IF (newlmt ) THEN
    1286           ierr = NF_PUT_VARA_DOUBLE (nid,id_FOCE,debut,epais
    1287      $        ,pctsrf_t(1,is_oce,k))
    1288           ierr = NF_PUT_VARA_DOUBLE (nid,id_FSIC,debut,epais
    1289      $        ,pctsrf_t(1,is_sic,k))
    1290           ierr = NF_PUT_VARA_DOUBLE (nid,id_FTER,debut,epais
    1291      $        ,pctsrf_t(1,is_ter,k))
    1292           ierr = NF_PUT_VARA_DOUBLE (nid,id_FLIC,debut,epais
    1293      $        ,pctsrf_t(1,is_lic,k))
    1294       ELSE
    1295           ierr = NF_PUT_VARA_DOUBLE (nid,id_NAT,debut,epais
    1296      $        ,phy_nat(1,k))
    1297       ENDIF
    1298 c
    1299       ierr = NF_PUT_VARA_DOUBLE (nid,id_SST,debut,epais,phy_sst(1,k))
    1300       ierr = NF_PUT_VARA_DOUBLE (nid,id_BILS,debut,epais,phy_bil(1,k))
    1301       ierr = NF_PUT_VARA_DOUBLE (nid,id_ALB,debut,epais,phy_alb(1,k))
    1302       ierr = NF_PUT_VARA_DOUBLE (nid,id_RUG,debut,epais,phy_rug(1,k))
    1303 #else
    1304       ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k))
    1305       IF (newlmt ) THEN
    1306           ierr = NF_PUT_VARA_REAL (nid,id_FOCE,debut,epais
    1307      $        ,pctsrf_t(1,is_oce,k))
    1308           ierr = NF_PUT_VARA_REAL (nid,id_FSIC,debut,epais
    1309      $        ,pctsrf_t(1,is_sic,k))
    1310           ierr = NF_PUT_VARA_REAL (nid,id_FTER,debut,epais
    1311      $        ,pctsrf_t(1,is_ter,k))
    1312           ierr = NF_PUT_VARA_REAL (nid,id_FLIC,debut,epais
    1313      $        ,pctsrf_t(1,is_lic,k))
    1314       ELSE
    1315           ierr = NF_PUT_VARA_REAL (nid,id_NAT,debut,epais
    1316      $        ,phy_nat(1,k))
    1317       ENDIF
    1318       ierr = NF_PUT_VARA_REAL (nid,id_SST,debut,epais,phy_sst(1,k))
    1319       ierr = NF_PUT_VARA_REAL (nid,id_BILS,debut,epais,phy_bil(1,k))
    1320       ierr = NF_PUT_VARA_REAL (nid,id_ALB,debut,epais,phy_alb(1,k))
    1321       ierr = NF_PUT_VARA_REAL (nid,id_RUG,debut,epais,phy_rug(1,k))
    1322 #endif
    1323 c
    1324       ENDDO
    1325 c
    1326       ierr = NF_CLOSE(nid)
    1327 c
    1328 #else
    1329       WRITE(lunout,*)
    1330      & 'limit_netcdf: Earth-specific routine, needs Earth physics'
    1331 #endif
    1332 ! of #ifdef CPP_EARTH
    1333       STOP
    1334       END
     306!-------------------------------------------------------------------------------
     307! Arguments:
     308  CHARACTER(LEN=*),  INTENT(IN)     :: fnam     ! NetCDF file name
     309  CHARACTER(LEN=3),  INTENT(IN)     :: mode     ! RUG, SIC, SST or ALB
     310  LOGICAL,           INTENT(IN)     :: ibar     ! interp on pressure levels
     311  INTEGER,           INTENT(IN)     :: ndays    ! current year number of days
     312  REAL,    POINTER,  DIMENSION(:,:) :: champo   ! output field = f(t)
     313  LOGICAL, OPTIONAL,                      INTENT(IN) :: flag
     314! flag=T means:   extrapolation (SST case)  or  old ice (SIC case)
     315  REAL,    OPTIONAL, DIMENSION(iim,jjp1), INTENT(IN) :: mask
     316!-------------------------------------------------------------------------------
     317! Local variables:
     318!--- NetCDF
     319  INTEGER :: ierr, ncid, varid                  ! NetCDF identifiers
     320  CHARACTER(LEN=30)               :: dnam       ! dimension name
     321  CHARACTER(LEN=80)               :: varname    ! NetCDF variable name
     322!--- dimensions
     323  INTEGER,           DIMENSION(4) :: dids       ! NetCDF dimensions identifiers
     324  REAL, ALLOCATABLE, DIMENSION(:) :: dlon_ini   ! initial longitudes vector
     325  REAL, ALLOCATABLE, DIMENSION(:) :: dlat_ini   ! initial latitudes  vector
     326  REAL, POINTER,     DIMENSION(:) :: dlon, dlat ! reordered lon/lat  vectors
     327!--- fields
     328  INTEGER :: imdep, jmdep, lmdep                ! dimensions of 'champ'
     329  REAL, ALLOCATABLE, DIMENSION(:) :: champ      ! wanted field on initial grid
     330  REAL, ALLOCATABLE, DIMENSION(:) :: yder, timeyear
     331  REAL,              DIMENSION(iim,jjp1) :: champint   ! interpolated field
     332  REAL, ALLOCATABLE, DIMENSION(:,:,:)    :: champtime
     333  REAL, ALLOCATABLE, DIMENSION(:,:,:)    :: champan
     334  REAL                              :: time, by
     335!--- input files
     336  CHARACTER(LEN=20)                 :: cal_in   ! calendar
     337  INTEGER                           :: ndays_in ! number of days
     338!--- misc
     339  INTEGER :: i, j, k, l                         ! loop counters
     340  REAL, ALLOCATABLE, DIMENSION(:,:) :: work     ! used for extrapolation
     341  CHARACTER(LEN=25)                 :: title    ! for messages
     342  LOGICAL                           :: extrp    ! flag for extrapolation
     343  REAL                              :: chmin, chmax
     344!-------------------------------------------------------------------------------
     345!---Variables depending on keyword 'mode' --------------------------------------
     346  NULLIFY(champo)
     347  SELECT CASE(mode)
     348    CASE('RUG'); varname='RUGOS';  title='Rugosite'
     349    CASE('SIC'); varname='sicbcs'; title='Sea-ice'
     350    CASE('SST'); varname='tosbcs'; title='SST'
     351    CASE('ALB'); varname='ALBEDO'; title='Albedo'
     352  END SELECT
     353  extrp=(flag.AND.mode=='SST')
     354
     355!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE ------------------------------
     356  ierr=NF90_OPEN(fnam,NF90_NOWRITE,ncid);                  CALL ncerr(ierr,fnam)
     357  ierr=NF90_INQ_VARID(ncid,varname,varid);                 CALL ncerr(ierr,fnam)
     358  ierr=NF90_INQUIRE_VARIABLE(ncid,varid,dimids=dids);      CALL ncerr(ierr,fnam)
     359
     360!--- Longitude
     361  ierr=NF90_INQUIRE_DIMENSION(ncid,dids(1),name=dnam,len=imdep)
     362  CALL ncerr(ierr,fnam); ALLOCATE(dlon_ini(imdep),dlon(imdep))
     363  ierr=NF90_INQ_VARID(ncid,dnam,varid);                    CALL ncerr(ierr,fnam)
     364  ierr=NF90_GET_VAR(ncid,varid,dlon_ini);                  CALL ncerr(ierr,fnam)
     365  WRITE(lunout,*) 'variable ', dnam, 'dimension ', imdep
     366
     367!--- Latitude
     368  ierr=NF90_INQUIRE_DIMENSION(ncid,dids(2),name=dnam,len=jmdep)
     369  CALL ncerr(ierr,fnam); ALLOCATE(dlat_ini(jmdep),dlat(jmdep))
     370  ierr=NF90_INQ_VARID(ncid,dnam,varid);                    CALL ncerr(ierr,fnam)
     371  ierr=NF90_GET_VAR(ncid,varid,dlat_ini);                  CALL ncerr(ierr,fnam)
     372  WRITE(lunout,*) 'variable ', dnam, 'dimension ', jmdep
     373
     374!--- Time (variable is not needed - it is rebuilt - but calendar is)
     375  ierr=NF90_INQUIRE_DIMENSION(ncid,dids(3),name=dnam,len=lmdep)
     376  CALL ncerr(ierr,fnam); ALLOCATE(timeyear(lmdep))
     377  ierr=NF90_INQ_VARID(ncid,dnam,varid);                    CALL ncerr(ierr,fnam)
     378  cal_in=' '
     379  ierr=NF90_GET_ATT(ncid,varid,'calendar',cal_in)
     380  IF(ierr/=NF90_NOERR) THEN
     381    SELECT CASE(mode)
     382      CASE('RUG','ALB'); cal_in='360d'
     383      CASE('SIC','SST'); cal_in='gregorian'
     384    END SELECT
     385    WRITE(lunout,*)'ATTENTION: variable ''time'' sans attribut ''calendrier'' d&
     386     &ans '//TRIM(fnam)//'. On choisit la valeur par defaut.'
     387  END IF
     388  WRITE(lunout,*) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', cal_in
     389
     390!--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION ----------------------
     391  !--- Determining input file number of days, depending on calendar
     392  ndays_in=year_len(anneeref,cal_in)
     393
     394!--- Time vector reconstruction (time vector from file is not trusted)
     395!--- If input records are not monthly, time sampling has to be constant !
     396  timeyear=mid_months(anneeref,cal_in,lmdep)
     397  IF(lmdep/=12) WRITE(lunout,'(a,i3,a)')'Note: les fichiers de '//TRIM(mode)   &
     398    //' ne comportent pas 12, mais ',lmdep,' renregistrements.'
     399
     400!--- GETTING THE FIELD AND INTERPOLATING IT ------------------------------------
     401  ALLOCATE(champ(imdep*jmdep),champtime(iim,jjp1,lmdep))
     402  IF(extrp) ALLOCATE(work(imdep,jmdep))
     403
     404  WRITE(lunout,*)
     405  WRITE(lunout,'(a,i3,a)')'LECTURE ET INTERPOLATION HORIZ. DE ',lmdep,' CHAMPS.'
     406  ierr=NF90_INQ_VARID(ncid,varname,varid);                 CALL ncerr(ierr,fnam)
     407  DO l=1,lmdep
     408    ierr=NF90_GET_VAR(ncid,varid,champ,(/1,1,l/),(/imdep,jmdep,1/))
     409    CALL ncerr(ierr,fnam)
     410    CALL conf_dat2d(title,imdep,jmdep,dlon_ini,dlat_ini,dlon,dlat,champ,ibar)
     411
     412    IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
     413
     414    IF(ibar.AND..NOT.(mode=='SIC'.AND.flag)) THEN
     415      IF(l==1) THEN
     416        WRITE(lunout,*)                                                        &
     417  '-------------------------------------------------------------------------'
     418        WRITE(lunout,*)                                                        &
     419  '$$$ Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$'
     420        WRITE(lunout,*)                                                        &
     421  '-------------------------------------------------------------------------'
     422      END IF
     423      IF(mode=='RUG') champ=LOG(champ)
     424      CALL inter_barxy(imdep, jmdep-1, dlon, dlat, champ, iim, jjm, rlonu,     &
     425                         rlatv, jjp1, champint)
     426      IF(mode=='RUG') THEN
     427        champint=EXP(champint)
     428        WHERE(NINT(mask)/=1) champint=0.001
     429      END IF
     430    ELSE
     431      SELECT CASE(mode)
     432        CASE('RUG');       CALL rugosite(imdep, jmdep, dlon, dlat, champ,      &
     433                                    iim, jjp1, rlonv, rlatu, champint, mask)
     434        CASE('SIC');       CALL sea_ice (imdep, jmdep, dlon, dlat, champ,      &
     435                                    iim, jjp1, rlonv, rlatu, champint)
     436        CASE('SST','ALB'); CALL grille_m(imdep, jmdep, dlon, dlat, champ,      &
     437                                    iim, jjp1, rlonv, rlatu, champint)
     438      END SELECT
     439    END IF
     440    champtime(:,:,l)=champint
     441  END DO
     442  ierr=NF90_CLOSE(ncid);                                   CALL ncerr(ierr,fnam)
     443
     444  DEALLOCATE(dlon_ini,dlat_ini,dlon,dlat,champ)
     445  IF(extrp) DEALLOCATE(work)
     446
     447!--- TIME INTERPOLATION --------------------------------------------------------
     448  WRITE(lunout,*)
     449  WRITE(lunout,*)'INTERPOLATION TEMPORELLE.'
     450  WRITE(lunout,"(2x,' Vecteur temps en entree: ',10f6.1)") timeyear
     451  WRITE(lunout,"(2x,' Vecteur temps en sortie de 0 a ',i3)") ndays
     452  ALLOCATE(yder(lmdep),champan(iip1,jjp1,ndays))
     453  DO j=1,jjp1
     454    DO i=1,iim
     455      CALL spline(timeyear,champtime(i,j,:),lmdep,1.e30,1.e30,yder)
     456      DO k=1,ndays
     457        time=FLOAT((k-1)*ndays_in)/FLOAT(ndays)
     458        CALL splint(timeyear,champtime(i,j,:),yder,lmdep,time,by)
     459        champan(i,j,k) = by
     460      END DO
     461    END DO
     462  END DO
     463  champan(iip1,:,:)=champan(1,:,:)
     464  DEALLOCATE(yder,champtime,timeyear)
     465
     466!--- Checking the result
     467  DO j=1,jjp1
     468    CALL minmax(iip1,champan(1,j,10),chmin,chmax)
     469    WRITE(lunout,*)' '//TRIM(title)//' au temps 10 ',chmin,chmax,j
     470  END DO
     471
     472!--- SPECIAL FILTER FOR SST: SST>271.38 ----------------------------------------
     473  IF(mode=='SST') THEN
     474    WRITE(lunout,*) 'Filtrage de la SST: SST >= 271.38'
     475    WHERE(champan<271.38) champan=271.38
     476  END IF
     477
     478!--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 ---------------------------------------
     479  IF(mode=='SIC') THEN
     480    WRITE(lunout,*) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
     481    champan(:,:,:)=champan(:,:,:)/100.
     482    champan(iip1,:,:)=champan(1,:,:)
     483    WHERE(champan>1.0) champan=1.0
     484    WHERE(champan<0.0) champan=0.0
     485  END IF
     486
     487write(*,*)'coin1'
     488!--- DYNAMICAL TO PHYSICAL GRID ------------------------------------------------
     489  ALLOCATE(champo(klon,ndays))
     490  DO k=1,ndays
     491    CALL gr_dyn_fi(1,iip1,jjp1,klon,champan(1,1,k),champo(1,k))
     492  END DO
     493write(*,*)'coin2'
     494  DEALLOCATE(champan)
     495write(*,*)'coin3'
     496END SUBROUTINE get_2Dfield
     497!
     498!-------------------------------------------------------------------------------
     499
     500
     501
     502!-------------------------------------------------------------------------------
     503!
     504FUNCTION year_len(y,cal_in)
     505!
     506!-------------------------------------------------------------------------------
     507  USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_year_len
     508  IMPLICIT NONE
     509!-------------------------------------------------------------------------------
     510! Arguments:
     511  INTEGER                       :: year_len
     512  INTEGER,           INTENT(IN) :: y
     513  CHARACTER(LEN=*),  INTENT(IN) :: cal_in
     514!-------------------------------------------------------------------------------
     515! Local variables:
     516  CHARACTER(LEN=20)             :: cal_out              ! calendar (for outputs)
     517!-------------------------------------------------------------------------------
     518!--- Getting the input calendar to reset at the end of the function
     519  CALL ioget_calendar(cal_out)
     520
     521!--- Unlocking calendar and setting it to wanted one
     522  CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
     523
     524!--- Getting the number of days in this year
     525  year_len=ioget_year_len(y)
     526
     527!--- Back to original calendar
     528  CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
     529
     530END FUNCTION year_len
     531!
     532!-------------------------------------------------------------------------------
     533
     534
     535!-------------------------------------------------------------------------------
     536!
     537FUNCTION mid_months(y,cal_in,nm)
     538!
     539!-------------------------------------------------------------------------------
     540  USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_mon_len
     541  IMPLICIT NONE
     542!-------------------------------------------------------------------------------
     543! Arguments:
     544  INTEGER,                INTENT(IN) :: y               ! year
     545  CHARACTER(LEN=*),       INTENT(IN) :: cal_in          ! calendar
     546  INTEGER,                INTENT(IN) :: nm              ! months/year number
     547  REAL,    DIMENSION(nm)             :: mid_months      ! mid-month times
     548!-------------------------------------------------------------------------------
     549! Local variables:
     550  CHARACTER(LEN=99)                  :: mess            ! error message
     551  CHARACTER(LEN=20)                  :: cal_out         ! calendar (for outputs)
     552  INTEGER, DIMENSION(nm)             :: mnth            ! months lengths (days)
     553  INTEGER                            :: m               ! months counter
     554  INTEGER                            :: nd              ! number of days
     555!-------------------------------------------------------------------------------
     556  nd=year_len(y,cal_in)
     557
     558  IF(nm==12) THEN
     559
     560  !--- Getting the input calendar to reset at the end of the function
     561    CALL ioget_calendar(cal_out)
     562
     563  !--- Unlocking calendar and setting it to wanted one
     564    CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
     565
     566  !--- Getting the length of each month
     567    DO m=1,nm; mnth(m)=ioget_mon_len(y,m); END DO
     568
     569  !--- Back to original calendar
     570    CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
     571
     572  ELSE IF(MODULO(nd,nm)/=0) THEN
     573    WRITE(mess,'(a,i3,a,i3,a)')'Unconsistent calendar: ',nd,' days/year, but ',&
     574      nm,' months/year. Months number should divide days number.'
     575    CALL abort_gcm('mid_months',TRIM(mess),1)
     576
     577  ELSE
     578    mnth=(/(m,m=1,nm,nd/nm)/)
     579  END IF
     580
     581!--- Mid-months times
     582  mid_months(1)=0.5*FLOAT(mnth(1))
     583  DO k=2,nm
     584    mid_months(k)=mid_months(k-1)+0.5*FLOAT(mnth(k-1)+mnth(k))
     585  END DO
     586
     587END FUNCTION mid_months
     588!
     589!-------------------------------------------------------------------------------
     590
     591
     592!-------------------------------------------------------------------------------
     593!
     594SUBROUTINE ncerr(ncres,fnam)
     595!
     596!-------------------------------------------------------------------------------
     597! Purpose: NetCDF errors handling.
     598!-------------------------------------------------------------------------------
     599  USE netcdf, ONLY : NF90_NOERR, NF90_STRERROR
     600  IMPLICIT NONE
     601!-------------------------------------------------------------------------------
     602! Arguments:
     603  INTEGER,          INTENT(IN) :: ncres
     604  CHARACTER(LEN=*), INTENT(IN) :: fnam
     605!-------------------------------------------------------------------------------
     606#include "iniprint.h"
     607  IF(ncres/=NF90_NOERR) THEN
     608    WRITE(lunout,*)'Problem with file '//TRIM(fnam)//' in routine limit_netcdf.'
     609    CALL abort_gcm('limit_netcdf',NF90_STRERROR(ncres),1)
     610  END IF
     611
     612END SUBROUTINE ncerr
     613!
     614!-------------------------------------------------------------------------------
     615
     616
     617END SUBROUTINE limit_netcdf
Note: See TracChangeset for help on using the changeset viewer.