My Project
 All Classes Files Functions Variables Macros
limit_netcdf.F90
Go to the documentation of this file.
1 !
2 ! $Id: limit_netcdf.F90 1508 2011-04-15 13:05:34Z jghattas $
3 !-------------------------------------------------------------------------------
4 !
5 SUBROUTINE 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 !-------------------------------------------------------------------------------
22  USE control_mod
23 #ifdef CPP_EARTH
24  USE dimphy
25  USE ioipsl, ONLY : ioget_year_len
26  USE phys_state_var_mod, ONLY : pctsrf
27  USE netcdf, ONLY : NF90_OPEN, NF90_CREATE, NF90_CLOSE, &
28  nf90_def_dim, nf90_def_var, nf90_put_var, nf90_put_att, &
29  nf90_noerr, nf90_nowrite, nf90_double, nf90_global, &
30  nf90_clobber, nf90_enddef, nf90_unlimited, nf90_float
31  USE inter_barxy_m, only: inter_barxy
32 #endif
33  IMPLICIT NONE
34 !-------------------------------------------------------------------------------
35 ! Arguments:
36 #include "dimensions.h"
37 #include "paramet.h"
38 #include "iniprint.h"
39  LOGICAL, INTENT(IN) :: interbar ! barycentric interpolation
40  LOGICAL, INTENT(IN) :: extrap ! SST extrapolation flag
41  LOGICAL, INTENT(IN) :: oldice ! old way ice computation
42  REAL, DIMENSION(iip1,jjp1), INTENT(IN) :: masque ! land mask
43 #ifndef CPP_EARTH
44  CALL abort_gcm('limit_netcdf','Earth-specific routine, needs Earth physics',1)
45 #else
46 !-------------------------------------------------------------------------------
47 ! Local variables:
48 #include "logic.h"
49 #include "comvert.h"
50 #include "comgeom2.h"
51 #include "comconst.h"
52 #include "indicesol.h"
53 
54 !--- INPUT NETCDF FILES NAMES --------------------------------------------------
55  CHARACTER(LEN=25) :: icefile, sstfile, dumstr
56  CHARACTER(LEN=25), PARAMETER :: famipsst='amipbc_sst_1x1.nc ', &
57  famipsic='amipbc_sic_1x1.nc ', &
58  fcpldsst='cpl_atm_sst.nc ', &
59  fcpldsic='cpl_atm_sic.nc ', &
60  fhistsst='histmth_sst.nc ', &
61  fhistsic='histmth_sic.nc ', &
62  frugo ='Rugos.nc ', &
63  falbe ='Albedo.nc '
64  CHARACTER(LEN=10) :: varname
65 !--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------
66  REAL, DIMENSION(klon) :: fi_ice, verif
67  REAL, DIMENSION(:,:), POINTER :: phy_rug=>null(), phy_ice=>null()
68  REAL, DIMENSION(:,:), POINTER :: phy_sst=>null(), phy_alb=>null()
69  REAL, DIMENSION(:,:), ALLOCATABLE :: phy_bil
70  REAL, DIMENSION(:,:,:), ALLOCATABLE :: pctsrf_t
71  INTEGER :: nbad
72 
73 !--- VARIABLES FOR OUTPUT FILE WRITING -----------------------------------------
74  INTEGER :: ierr, nid, ndim, ntim, k
75  INTEGER, DIMENSION(2) :: dims
76  INTEGER :: id_tim, id_sst, id_bils, id_rug, id_alb
77  INTEGER :: id_foce, id_fsic, id_fter, id_flic
78  INTEGER :: nf90_format
79  INTEGER :: ndays !--- Depending on the output calendar
80 
81 !--- INITIALIZATIONS -----------------------------------------------------------
82 #ifdef NC_DOUBLE
83  nf90_format=nf90_double
84 #else
85  nf90_format=nf90_float
86 #endif
87 
88  pi = 4.*atan(1.)
89  rad = 6371229.
90  daysec= 86400.
91  omeg = 2.*pi/daysec
92  g = 9.8
93  kappa = 0.2857143
94  cpp = 1004.70885
95  dtvr = daysec/REAL(day_step)
96  CALL inigeom
97 
98 !--- Beware: anneeref (from gcm.def) is used to determine output time sampling
99  ndays=ioget_year_len(anneeref)
100 
101 !--- RUGOSITY TREATMENT --------------------------------------------------------
102  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la rugosite'
103  varname='RUGOS'
104  CALL get_2dfield(frugo,varname,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:))
105 
106 !--- OCEAN TREATMENT -----------------------------------------------------------
107  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la glace oceanique'
108 
109 ! Input SIC file selection
110 ! Open file only to test if available
111  IF ( nf90_open(trim(famipsic),nf90_nowrite,nid)==nf90_noerr ) THEN
112  icefile=trim(famipsic)
113  varname='sicbcs'
114  ELSE IF( nf90_open(trim(fcpldsic),nf90_nowrite,nid)==nf90_noerr ) THEN
115  icefile=trim(fcpldsic)
116  varname='SIICECOV'
117  ELSE IF ( nf90_open(trim(fhistsic),nf90_nowrite,nid)==nf90_noerr ) THEN
118  icefile=trim(fhistsic)
119  varname='pourc_sic'
120  ELSE
121  WRITE(lunout,*) 'ERROR! No sea-ice input file was found.'
122  WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsic),', ',trim(fcpldsic),', ',trim(fhistsic)
123  CALL abort_gcm('limit_netcdf','No sea-ice file was found',1)
124  END IF
125  ierr=nf90_close(nid)
126  IF (prt_level>=0) WRITE(lunout,*)'Pour la glace de mer a ete choisi le fichier '//trim(icefile)
127 
128  CALL get_2dfield(icefile,varname, 'SIC',interbar,ndays,phy_ice,flag=oldice)
129 
130  ALLOCATE(pctsrf_t(klon,nbsrf,ndays))
131  DO k=1,ndays
132  fi_ice=phy_ice(:,k)
133  WHERE(fi_ice>=1.0 ) fi_ice=1.0
134  WHERE(fi_ice<epsfra) fi_ice=0.0
135  pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter) ! land soil
136  pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic) ! land ice
137  IF (icefile==trim(fcpldsic)) THEN ! SIC=pICE*(1-LIC-TER)
138  pctsrf_t(:,is_sic,k)=fi_ice(:)*(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter))
139  ELSE IF (icefile==trim(fhistsic)) THEN ! SIC=pICE
140  pctsrf_t(:,is_sic,k)=fi_ice(:)
141  ELSE ! icefile==famipsic ! SIC=pICE-LIC
142  pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k)
143  END IF
144  WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0.
145  WHERE(1.0-zmasq<epsfra)
146  pctsrf_t(:,is_sic,k)=0.0
147  pctsrf_t(:,is_oce,k)=0.0
148  ELSEWHERE
149  WHERE(pctsrf_t(:,is_sic,k)>=1.0-zmasq)
150  pctsrf_t(:,is_sic,k)=1.0-zmasq
151  pctsrf_t(:,is_oce,k)=0.0
152  ELSEWHERE
153  pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)
154  WHERE(pctsrf_t(:,is_oce,k)<epsfra)
155  pctsrf_t(:,is_oce,k)=0.0
156  pctsrf_t(:,is_sic,k)=1.0-zmasq
157  END WHERE
158  END WHERE
159  END WHERE
160  nbad=count(pctsrf_t(:,is_oce,k)<0.0)
161  IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
162  nbad=count(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>epsfra)
163  IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad
164  END DO
165  DEALLOCATE(phy_ice)
166 
167 !--- SST TREATMENT -------------------------------------------------------------
168  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la sst'
169 
170 ! Input SST file selection
171 ! Open file only to test if available
172  IF ( nf90_open(trim(famipsst),nf90_nowrite,nid)==nf90_noerr ) THEN
173  sstfile=trim(famipsst)
174  varname='tosbcs'
175  ELSE IF ( nf90_open(trim(fcpldsst),nf90_nowrite,nid)==nf90_noerr ) THEN
176  sstfile=trim(fcpldsst)
177  varname='SISUTESW'
178  ELSE IF ( nf90_open(trim(fhistsst),nf90_nowrite,nid)==nf90_noerr ) THEN
179  sstfile=trim(fhistsst)
180  varname='tsol_oce'
181  ELSE
182  WRITE(lunout,*) 'ERROR! No sst input file was found.'
183  WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsst),trim(fcpldsst),trim(fhistsst)
184  CALL abort_gcm('limit_netcdf','No sst file was found',1)
185  END IF
186  ierr=nf90_close(nid)
187  IF (prt_level>=0) WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//trim(sstfile)
188 
189  CALL get_2dfield(sstfile,varname,'SST',interbar,ndays,phy_sst,flag=extrap)
190 
191 !--- ALBEDO TREATMENT ----------------------------------------------------------
192  IF (prt_level>1) WRITE(lunout,*) 'Traitement de l albedo'
193  varname='ALBEDO'
194  CALL get_2dfield(falbe,varname,'ALB',interbar,ndays,phy_alb)
195 
196 !--- REFERENCE GROUND HEAT FLUX TREATMENT --------------------------------------
197  ALLOCATE(phy_bil(klon,ndays)); phy_bil=0.0
198 
199 !--- OUTPUT FILE WRITING -------------------------------------------------------
200  IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : debut'
201 
202  !--- File creation
203  ierr=nf90_create("limit.nc",nf90_clobber,nid)
204  ierr=nf90_put_att(nid,nf90_global,"title","Fichier conditions aux limites")
205 
206  !--- Dimensions creation
207  ierr=nf90_def_dim(nid,"points_physiques",klon,ndim)
208  ierr=nf90_def_dim(nid,"time",nf90_unlimited,ntim)
209 
210  dims=(/ndim,ntim/)
211 
212  !--- Variables creation
213  ierr=nf90_def_var(nid,"TEMPS",nf90_format,(/ntim/),id_tim)
214  ierr=nf90_def_var(nid,"FOCE", nf90_format,dims,id_foce)
215  ierr=nf90_def_var(nid,"FSIC", nf90_format,dims,id_fsic)
216  ierr=nf90_def_var(nid,"FTER", nf90_format,dims,id_fter)
217  ierr=nf90_def_var(nid,"FLIC", nf90_format,dims,id_flic)
218  ierr=nf90_def_var(nid,"SST", nf90_format,dims,id_sst)
219  ierr=nf90_def_var(nid,"BILS", nf90_format,dims,id_bils)
220  ierr=nf90_def_var(nid,"ALB", nf90_format,dims,id_alb)
221  ierr=nf90_def_var(nid,"RUG", nf90_format,dims,id_rug)
222 
223  !--- Attributes creation
224  ierr=nf90_put_att(nid,id_tim, "title","Jour dans l annee")
225  ierr=nf90_put_att(nid,id_foce,"title","Fraction ocean")
226  ierr=nf90_put_att(nid,id_fsic,"title","Fraction glace de mer")
227  ierr=nf90_put_att(nid,id_fter,"title","Fraction terre")
228  ierr=nf90_put_att(nid,id_flic,"title","Fraction land ice")
229  ierr=nf90_put_att(nid,id_sst ,"title","Temperature superficielle de la mer")
230  ierr=nf90_put_att(nid,id_bils,"title","Reference flux de chaleur au sol")
231  ierr=nf90_put_att(nid,id_alb, "title","Albedo a la surface")
232  ierr=nf90_put_att(nid,id_rug, "title","Rugosite")
233 
234  ierr=nf90_enddef(nid)
235 
236  !--- Variables saving
237  ierr=nf90_put_var(nid,id_tim,(/(REAL(k),k=1,ndays)/))
238  ierr=nf90_put_var(nid,id_foce,pctsrf_t(:,is_oce,:),(/1,1/),(/klon,ndays/))
239  ierr=nf90_put_var(nid,id_fsic,pctsrf_t(:,is_sic,:),(/1,1/),(/klon,ndays/))
240  ierr=nf90_put_var(nid,id_fter,pctsrf_t(:,is_ter,:),(/1,1/),(/klon,ndays/))
241  ierr=nf90_put_var(nid,id_flic,pctsrf_t(:,is_lic,:),(/1,1/),(/klon,ndays/))
242  ierr=nf90_put_var(nid,id_sst ,phy_sst(:,:),(/1,1/),(/klon,ndays/))
243  ierr=nf90_put_var(nid,id_bils,phy_bil(:,:),(/1,1/),(/klon,ndays/))
244  ierr=nf90_put_var(nid,id_alb ,phy_alb(:,:),(/1,1/),(/klon,ndays/))
245  ierr=nf90_put_var(nid,id_rug ,phy_rug(:,:),(/1,1/),(/klon,ndays/))
246 
247  ierr=nf90_close(nid)
248 
249  IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : fin'
250 
251  DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug)
252 
253 
254 !===============================================================================
255 !
256  CONTAINS
257 !
258 !===============================================================================
259 
260 
261 !-------------------------------------------------------------------------------
262 !
263 SUBROUTINE get_2dfield(fnam, varname, mode, ibar, ndays, champo, flag, mask)
264 !
265 !-----------------------------------------------------------------------------
266 ! Comments:
267 ! There are two assumptions concerning the NetCDF files, that are satisfied
268 ! with files that are conforming NC convention:
269 ! 1) The last dimension of the variables used is the time record.
270 ! 2) Dimensional variables have the same names as corresponding dimensions.
271 !-----------------------------------------------------------------------------
272  USE netcdf, ONLY: NF90_OPEN, NF90_INQ_VARID, NF90_INQUIRE_VARIABLE, &
273  nf90_close, nf90_inq_dimid, nf90_inquire_dimension, nf90_get_var, &
274  nf90_get_att
275  USE dimphy, ONLY : klon
276  USE phys_state_var_mod, ONLY : pctsrf
277  USE control_mod
278  use pchsp_95_m, only: pchsp_95
279  use pchfe_95_m, only: pchfe_95
280  use arth_m, only: arth
281 
282  IMPLICIT NONE
283 #include "dimensions.h"
284 #include "paramet.h"
285 #include "comgeom2.h"
286 #include "indicesol.h"
287 #include "iniprint.h"
288 !-----------------------------------------------------------------------------
289 ! Arguments:
290  CHARACTER(LEN=*), INTENT(IN) :: fnam ! NetCDF file name
291  CHARACTER(LEN=10), INTENT(IN) :: varname ! NetCDF variable name
292  CHARACTER(LEN=3), INTENT(IN) :: mode ! RUG, SIC, SST or ALB
293  LOGICAL, INTENT(IN) :: ibar ! interp on pressure levels
294  INTEGER, INTENT(IN) :: ndays ! current year number of days
295  REAL, POINTER, DIMENSION(:, :) :: champo ! output field = f(t)
296  LOGICAL, OPTIONAL, INTENT(IN) :: flag ! extrapol. (SST) old ice (SIC)
297  REAL, OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask
298 !------------------------------------------------------------------------------
299 ! Local variables:
300 !--- NetCDF
301  INTEGER :: ncid, varid ! NetCDF identifiers
302  CHARACTER(LEN=30) :: dnam ! dimension name
303 !--- dimensions
304  INTEGER, DIMENSION(4) :: dids ! NetCDF dimensions identifiers
305  REAL, ALLOCATABLE, DIMENSION(:) :: dlon_ini ! initial longitudes vector
306  REAL, ALLOCATABLE, DIMENSION(:) :: dlat_ini ! initial latitudes vector
307  REAL, POINTER, DIMENSION(:) :: dlon, dlat ! reordered lon/lat vectors
308 !--- fields
309  INTEGER :: imdep, jmdep, lmdep ! dimensions of 'champ'
310  REAL, ALLOCATABLE, DIMENSION(:, :) :: champ ! wanted field on initial grid
311  REAL, ALLOCATABLE, DIMENSION(:) :: yder, timeyear
312  REAL, DIMENSION(iim, jjp1) :: champint ! interpolated field
313  REAL, ALLOCATABLE, DIMENSION(:, :, :) :: champtime
314  REAL, ALLOCATABLE, DIMENSION(:, :, :) :: champan
315 !--- input files
316  CHARACTER(LEN=20) :: cal_in ! calendar
317  CHARACTER(LEN=20) :: unit_sic ! attribute unit in sea-ice file
318  INTEGER :: ndays_in ! number of days
319 !--- misc
320  INTEGER :: i, j, k, l ! loop counters
321  REAL, ALLOCATABLE, DIMENSION(:, :) :: work ! used for extrapolation
322  CHARACTER(LEN=25) :: title ! for messages
323  LOGICAL :: extrp ! flag for extrapolation
324  LOGICAL :: oldice ! flag for old way ice computation
325  REAL :: chmin, chmax
326  INTEGER ierr
327  integer n_extrap ! number of extrapolated points
328  logical skip
329 
330 !------------------------------------------------------------------------------
331 !---Variables depending on keyword 'mode' -------------------------------------
332  nullify(champo)
333 
334  SELECT CASE(mode)
335  CASE('RUG'); title='Rugosite'
336  CASE('SIC'); title='Sea-ice'
337  CASE('SST'); title='SST'
338  CASE('ALB'); title='Albedo'
339  END SELECT
340 
341 
342  extrp=.false.
343  oldice=.false.
344  IF ( present(flag) ) THEN
345  IF ( flag .AND. mode=='SST' ) extrp=.true.
346  IF ( flag .AND. mode=='SIC' ) oldice=.true.
347  END IF
348 
349 !--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
350  IF (prt_level>5) WRITE(lunout,*) ' Now reading file : ',fnam
351  ierr=nf90_open(fnam, nf90_nowrite, ncid); CALL ncerr(ierr, fnam)
352  ierr=nf90_inq_varid(ncid, trim(varname), varid); CALL ncerr(ierr, fnam)
353  ierr=nf90_inquire_variable(ncid, varid, dimids=dids); CALL ncerr(ierr, fnam)
354 
355 !--- Read unit for sea-ice variable only
356  IF (mode=='SIC') THEN
357  ierr=nf90_get_att(ncid, varid, 'units', unit_sic)
358  IF(ierr/=nf90_noerr) THEN
359  IF (prt_level>5) WRITE(lunout,*) 'No unit was given in sea-ice file. Take percentage as default value'
360  unit_sic='X'
361  ELSE
362  IF (prt_level>5) WRITE(lunout,*) ' Sea-ice cover has unit=',unit_sic
363  END IF
364  END IF
365 
366 !--- Longitude
367  ierr=nf90_inquire_dimension(ncid, dids(1), name=dnam, len=imdep)
368  CALL ncerr(ierr, fnam); ALLOCATE(dlon_ini(imdep), dlon(imdep))
369  ierr=nf90_inq_varid(ncid, dnam, varid); CALL ncerr(ierr, fnam)
370  ierr=nf90_get_var(ncid, varid, dlon_ini); CALL ncerr(ierr, fnam)
371  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep
372 
373 !--- Latitude
374  ierr=nf90_inquire_dimension(ncid, dids(2), name=dnam, len=jmdep)
375  CALL ncerr(ierr, fnam); ALLOCATE(dlat_ini(jmdep), dlat(jmdep))
376  ierr=nf90_inq_varid(ncid, dnam, varid); CALL ncerr(ierr, fnam)
377  ierr=nf90_get_var(ncid, varid, dlat_ini); CALL ncerr(ierr, fnam)
378  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep
379 
380 !--- Time (variable is not needed - it is rebuilt - but calendar is)
381  ierr=nf90_inquire_dimension(ncid, dids(3), name=dnam, len=lmdep)
382  CALL ncerr(ierr, fnam); ALLOCATE(timeyear(lmdep))
383  ierr=nf90_inq_varid(ncid, dnam, varid); CALL ncerr(ierr, fnam)
384  cal_in=' '
385  ierr=nf90_get_att(ncid, varid, 'calendar', cal_in)
386  IF(ierr/=nf90_noerr) THEN
387  SELECT CASE(mode)
388  CASE('RUG', 'ALB'); cal_in='360d'
389  CASE('SIC', 'SST'); cal_in='gregorian'
390  END SELECT
391  IF (prt_level>5) WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' &
392  // 'dans '//trim(fnam)//'. On choisit la valeur par defaut.'
393  END IF
394  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', &
395  cal_in
396 
397 
398 !--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION --------------------
399  !--- Determining input file number of days, depending on calendar
400  ndays_in=year_len(anneeref, cal_in)
401 
402 !--- Time vector reconstruction (time vector from file is not trusted)
403 !--- If input records are not monthly, time sampling has to be constant !
404  timeyear=mid_months(anneeref, cal_in, lmdep)
405  IF (lmdep /= 12) WRITE(lunout,*) 'Note : les fichiers de ', trim(mode), &
406  ' ne comportent pas 12, mais ', lmdep, ' enregistrements.'
407 
408 !--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
409  ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep))
410  IF(extrp) ALLOCATE(work(imdep, jmdep))
411 
412  IF (prt_level>5) WRITE(lunout, *)
413  IF (prt_level>5) WRITE(lunout,*)'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, ' CHAMPS.'
414  ierr=nf90_inq_varid(ncid, varname, varid); CALL ncerr(ierr, fnam)
415  DO l=1, lmdep
416  ierr=nf90_get_var(ncid, varid, champ, (/1, 1, l/), (/imdep, jmdep, 1/))
417  CALL ncerr(ierr, fnam)
418  CALL conf_dat2d(title, imdep, jmdep, dlon_ini, dlat_ini, dlon, dlat, &
419  champ, ibar)
420 
421  IF (extrp) CALL extrapol(champ, imdep, jmdep, 999999., .true., .true., 2, &
422  work)
423 
424  IF(ibar .AND. .NOT.oldice) THEN
425  IF(l==1 .AND. prt_level>5) THEN
426  WRITE(lunout, *) '-------------------------------------------------------------------------'
427  WRITE(lunout, *) 'Utilisation de l''interpolation barycentrique pour ',trim(title),' $$$'
428  WRITE(lunout, *) '-------------------------------------------------------------------------'
429  END IF
430  IF(mode=='RUG') champ=log(champ)
431  CALL inter_barxy(dlon, dlat(:jmdep-1), champ, rlonu(:iim), &
432  rlatv, champint)
433  IF(mode=='RUG') THEN
434  champint=exp(champint)
435  WHERE(nint(mask)/=1) champint=0.001
436  END IF
437  ELSE
438  SELECT CASE(mode)
439  CASE('RUG'); CALL rugosite(imdep, jmdep, dlon, dlat, champ, &
440  iim, jjp1, rlonv, rlatu, champint, mask)
441  CASE('SIC'); CALL sea_ice(imdep, jmdep, dlon, dlat, champ, &
442  iim, jjp1, rlonv, rlatu, champint)
443  CASE('SST', 'ALB'); CALL grille_m(imdep, jmdep, dlon, dlat, champ, &
444  iim, jjp1, rlonv, rlatu, champint)
445  END SELECT
446  END IF
447  champtime(:, :, l)=champint
448  END DO
449  ierr=nf90_close(ncid); CALL ncerr(ierr, fnam)
450 
451  DEALLOCATE(dlon_ini, dlat_ini, dlon, dlat, champ)
452  IF(extrp) DEALLOCATE(work)
453 
454 !--- TIME INTERPOLATION ------------------------------------------------------
455  IF (prt_level>5) THEN
456  WRITE(lunout, *)
457  WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
458  WRITE(lunout, *)' Vecteur temps en entree: ', timeyear
459  WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays
460  END IF
461 
462  ALLOCATE(yder(lmdep), champan(iip1, jjp1, ndays))
463  skip = .false.
464  n_extrap = 0
465  DO j=1, jjp1
466  DO i=1, iim
467  yder = pchsp_95(timeyear, champtime(i, j, :), ibeg=2, iend=2, &
468  vc_beg=0., vc_end=0.)
469  CALL pchfe_95(timeyear, champtime(i, j, :), yder, skip, &
470  arth(0., real(ndays_in) / ndays, ndays), champan(i, j, :), ierr)
471  if (ierr < 0) stop 1
472  n_extrap = n_extrap + ierr
473  END DO
474  END DO
475  if (n_extrap /= 0) then
476  WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
477  end if
478  champan(iip1, :, :)=champan(1, :, :)
479  DEALLOCATE(yder, champtime, timeyear)
480 
481 !--- Checking the result
482  DO j=1, jjp1
483  CALL minmax(iip1, champan(1, j, 10), chmin, chmax)
484  IF (prt_level>5) WRITE(lunout, *)' ',trim(title),' au temps 10 ', chmin, chmax, j
485  END DO
486 
487 !--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
488  IF(mode=='SST') THEN
489  IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38'
490  WHERE(champan<271.38) champan=271.38
491  END IF
492 
493 !--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
494  IF(mode=='SIC') THEN
495  IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
496 
497  IF (unit_sic=='1') THEN
498  ! Nothing to be done for sea-ice field is already in fraction of 1
499  ! This is the case for sea-ice in file cpl_atm_sic.nc
500  IF (prt_level>5) WRITE(lunout,*) 'Sea-ice field already in fraction of 1'
501  ELSE
502  ! Convert sea ice from percentage to fraction of 1
503  IF (prt_level>5) WRITE(lunout,*) 'Transformt sea-ice field from percentage to fraction of 1.'
504  champan(:, :, :)=champan(:, :, :)/100.
505  END IF
506 
507  champan(iip1, :, :)=champan(1, :, :)
508  WHERE(champan>1.0) champan=1.0
509  WHERE(champan<0.0) champan=0.0
510  END IF
511 
512 !--- DYNAMICAL TO PHYSICAL GRID ----------------------------------------------
513  ALLOCATE(champo(klon, ndays))
514  DO k=1, ndays
515  CALL gr_dyn_fi(1, iip1, jjp1, klon, champan(1, 1, k), champo(1, k))
516  END DO
517  DEALLOCATE(champan)
518 
519 END SUBROUTINE get_2dfield
520 !
521 !-------------------------------------------------------------------------------
522 
523 
524 
525 !-------------------------------------------------------------------------------
526 !
527 FUNCTION year_len(y,cal_in)
528 !
529 !-------------------------------------------------------------------------------
530  USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_year_len
531  IMPLICIT NONE
532 !-------------------------------------------------------------------------------
533 ! Arguments:
534  INTEGER :: year_len
535  INTEGER, INTENT(IN) :: y
536  CHARACTER(LEN=*), INTENT(IN) :: cal_in
537 !-------------------------------------------------------------------------------
538 ! Local variables:
539  CHARACTER(LEN=20) :: cal_out ! calendar (for outputs)
540 !-------------------------------------------------------------------------------
541 !--- Getting the input calendar to reset at the end of the function
542  CALL ioget_calendar(cal_out)
543 
544 !--- Unlocking calendar and setting it to wanted one
545  CALL lock_calendar(.false.); CALL ioconf_calendar(trim(cal_in))
546 
547 !--- Getting the number of days in this year
548  year_len=ioget_year_len(y)
549 
550 !--- Back to original calendar
551  CALL lock_calendar(.false.); CALL ioconf_calendar(trim(cal_out))
552 
553 END FUNCTION year_len
554 !
555 !-------------------------------------------------------------------------------
556 
557 
558 !-------------------------------------------------------------------------------
559 !
560 FUNCTION mid_months(y,cal_in,nm)
561 !
562 !-------------------------------------------------------------------------------
563  USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_mon_len
564  IMPLICIT NONE
565 !-------------------------------------------------------------------------------
566 ! Arguments:
567  INTEGER, INTENT(IN) :: y ! year
568  CHARACTER(LEN=*), INTENT(IN) :: cal_in ! calendar
569  INTEGER, INTENT(IN) :: nm ! months/year number
570  REAL, DIMENSION(nm) :: mid_months ! mid-month times
571 !-------------------------------------------------------------------------------
572 ! Local variables:
573  CHARACTER(LEN=99) :: mess ! error message
574  CHARACTER(LEN=20) :: cal_out ! calendar (for outputs)
575  INTEGER, DIMENSION(nm) :: mnth ! months lengths (days)
576  INTEGER :: m ! months counter
577  INTEGER :: nd ! number of days
578 !-------------------------------------------------------------------------------
579  nd=year_len(y,cal_in)
580 
581  IF(nm==12) THEN
582 
583  !--- Getting the input calendar to reset at the end of the function
584  CALL ioget_calendar(cal_out)
585 
586  !--- Unlocking calendar and setting it to wanted one
587  CALL lock_calendar(.false.); CALL ioconf_calendar(trim(cal_in))
588 
589  !--- Getting the length of each month
590  DO m=1,nm; mnth(m)=ioget_mon_len(y,m); END DO
591 
592  !--- Back to original calendar
593  CALL lock_calendar(.false.); CALL ioconf_calendar(trim(cal_out))
594 
595  ELSE IF(modulo(nd,nm)/=0) THEN
596  WRITE(mess,'(a,i3,a,i3,a)')'Unconsistent calendar: ',nd,' days/year, but ',&
597  nm,' months/year. Months number should divide days number.'
598  CALL abort_gcm('mid_months',trim(mess),1)
599 
600  ELSE
601  mnth=(/(m,m=1,nm,nd/nm)/)
602  END IF
603 
604 !--- Mid-months times
605  mid_months(1)=0.5*REAL(mnth(1))
606  DO k=2,nm
607  mid_months(k)=mid_months(k-1)+0.5*REAL(mnth(k-1)+mnth(k))
608  END DO
609 
610 END FUNCTION mid_months
611 !
612 !-------------------------------------------------------------------------------
613 
614 
615 !-------------------------------------------------------------------------------
616 !
617 SUBROUTINE ncerr(ncres,fnam)
618 !
619 !-------------------------------------------------------------------------------
620 ! Purpose: NetCDF errors handling.
621 !-------------------------------------------------------------------------------
622  USE netcdf, ONLY : NF90_NOERR, NF90_STRERROR
623  IMPLICIT NONE
624 !-------------------------------------------------------------------------------
625 ! Arguments:
626  INTEGER, INTENT(IN) :: ncres
627  CHARACTER(LEN=*), INTENT(IN) :: fnam
628 !-------------------------------------------------------------------------------
629 #include "iniprint.h"
630  IF(ncres/=nf90_noerr) THEN
631  WRITE(lunout,*)'Problem with file '//trim(fnam)//' in routine limit_netcdf.'
632  CALL abort_gcm('limit_netcdf',nf90_strerror(ncres),1)
633  END IF
634 
635 END SUBROUTINE ncerr
636 !
637 !-------------------------------------------------------------------------------
638 
639 #endif
640 ! of #ifdef CPP_EARTH
641 
642 END SUBROUTINE limit_netcdf