Context Navigation

limit_netcdf.F90 @ 1319

Last change on this file since 1319 was 1319, checked in by Laurent Fairhead, 14 years ago

Modifications to the start and limit creation routines to account for different

calendars

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

surface fractions read in the limit file

Force readaerosol.F90 to read in aerosols file with 12 timesteps

Modifications aux routines de créations des fichiers start et limit pour prendre

en compte différents calendriers

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

compatible avec les fractions de surface lu dans le fichier limit

Forcer readaerosol à ne lire que des fichiers à 12 pas de temps

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 25.6 KB

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

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

Download in other formats:

Original Format