Context Navigation

limit_netcdf.F90 @ 2154

Last change on this file since 2154 was 2154, checked in by acozic, 10 years ago
To create limit.nc add the possibility to read sst and sic from erai files (wait as sst.nc and ci.nc by the model)
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory
File size: 27.5 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format