source: LMDZ4/trunk/libf/dyn3d/limit_netcdf.F90 @ 2061

Last change on this file since 2061 was 1425, checked in by lguez, 14 years ago

Replaced Numerical Recipes procedures for spline interpolation (not in
the public domain) by procedures from the Pchip package in the Slatec
library. This only affects the program "ce0l", not the program
"gcm". Tested on Brodie SX8 with "-debug" and "-prod", "-parallel
none" and "-parallel mpi". "start.nc" and "limit.nc" are
changed. "startphy.nc" is not changed. The relative change is of order
1e-7 or less. The revision makes the program faster (tested on Brodie
with "-prod -d 144x142x39", CPU time is 38 s, instead of 54
s). Procedures from Slatec are untouched, except for
"i1mach.F". Created procedures "pchfe_95" and "pchsp_95" which are
wrappers for "pchfe" and "pchsp" from Slatec. "pchfe_95" and
"pchsp_95" have a safer and simpler interface.

Replaced "make" by "sxgmake" in "arch-SX8_BRODIE.fcm". Added files for
compilation by FCM with "g95".

In "arch-linux-32bit.fcm", replaced "pgf90" by "pgf95". There was no
difference between "dev" and "debug" so added "-O1" to "dev". Added
debugging options. Removed "-Wl,-Bstatic
-L/usr/lib/gcc-lib/i386-linux/2.95.2", which usually produces an error
at link-time.

Bash is now ubiquitous while KornShell? is not so use Bash instead of
KornShell? in FCM.

Replaced some statements "write(6,*)" by "write(lunout,*)". Replaced
"stop" by "stop 1" in the case where "abort_gcm" is called with "ierr
/= 0". Removed "stop" statements at the end of procedures
"limit_netcdf" and main program "ce0l" (why not let the program end
normally?).

Made some arrays automatic instead of allocatable in "start_inter_3d".

Zeroed "wake_pe", "fm_therm", "entr_therm" and "detr_therm" in
"dyn3dpar/etat0_netcdf.F90". The parallel and sequential results of
"ce0l" are thus identical.

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