source: LMDZ5/branches/testing/libf/phylmd/limit_netcdf.F90 @ 1856

Last change on this file since 1856 was 1795, checked in by Ehouarn Millour, 11 years ago

Version testing basee sur la r1794


Testing release based on r1794

File size: 26.4 KB
Line 
1!
2! $Id: limit_netcdf.F90 1508 2011-04-15 13:05:34Z jghattas $
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  USE indice_sol_mod
24#ifdef CPP_EARTH
25  USE dimphy
26  USE ioipsl,             ONLY : ioget_year_len
27  USE phys_state_var_mod, ONLY : pctsrf
28  USE netcdf,             ONLY : NF90_OPEN,    NF90_CREATE,  NF90_CLOSE,       &
29                   NF90_DEF_DIM, NF90_DEF_VAR, NF90_PUT_VAR, NF90_PUT_ATT,     &
30                   NF90_NOERR,   NF90_NOWRITE, NF90_DOUBLE,  NF90_GLOBAL,      &
31                   NF90_CLOBBER, NF90_ENDDEF,  NF90_UNLIMITED, NF90_FLOAT
32  USE inter_barxy_m, only: inter_barxy
33#endif
34  IMPLICIT NONE
35!-------------------------------------------------------------------------------
36! Arguments:
37#include "dimensions.h"
38#include "paramet.h"
39#include "iniprint.h"
40  LOGICAL,                    INTENT(IN) :: interbar ! barycentric interpolation
41  LOGICAL,                    INTENT(IN) :: extrap   ! SST extrapolation flag
42  LOGICAL,                    INTENT(IN) :: oldice   ! old way ice computation
43  REAL, DIMENSION(iip1,jjp1), INTENT(IN) :: masque   ! land mask
44#ifndef CPP_EARTH
45  CALL abort_gcm('limit_netcdf','Earth-specific routine, needs Earth physics',1)
46#else
47!-------------------------------------------------------------------------------
48! Local variables:
49#include "logic.h"
50#include "comvert.h"
51#include "comgeom2.h"
52#include "comconst.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!
263SUBROUTINE 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  USE indice_sol_mod
282
283  IMPLICIT NONE
284#include "dimensions.h"
285#include "paramet.h"
286#include "comgeom2.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
519END SUBROUTINE get_2Dfield
520!
521!-------------------------------------------------------------------------------
522
523
524
525!-------------------------------------------------------------------------------
526!
527FUNCTION 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
553END FUNCTION year_len
554!
555!-------------------------------------------------------------------------------
556
557
558!-------------------------------------------------------------------------------
559!
560FUNCTION 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
610END FUNCTION mid_months
611!
612!-------------------------------------------------------------------------------
613
614
615!-------------------------------------------------------------------------------
616!
617SUBROUTINE 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
635END SUBROUTINE ncerr
636!
637!-------------------------------------------------------------------------------
638
639#endif
640! of #ifdef CPP_EARTH
641
642END SUBROUTINE limit_netcdf
Note: See TracBrowser for help on using the repository browser.