Context Navigation

limit_netcdf.F90 @ 1354

Last change on this file since 1354 was 1328, checked in by Laurent Fairhead, 15 years ago

Continent/ocean mask was transformed incorrectly from o2a.nc file
Change of unit for seaice depending on whether the initial data is from
climatology or model results

Erreur sur la transformation vers la grille physique du masque lu dans o2a.nc
Pas de division par 100 de la glace océanique si on lit des champs du modèle
couplé

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

File size: 25.5 KB

Line
1	!
2	! $Id: limit_netcdf.F90 1328 2010-03-18 13:26:23Z musat $
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	USE inter_barxy_m, only: inter_barxy
31	#endif
32	IMPLICIT NONE
33	!-------------------------------------------------------------------------------
34	! Arguments:
35	#include "dimensions.h"
36	#include "paramet.h"
37	#include "iniprint.h"
38	LOGICAL, INTENT(IN) :: interbar ! barycentric interpolation
39	LOGICAL, INTENT(IN) :: extrap ! SST extrapolation flag
40	LOGICAL, INTENT(IN) :: oldice ! old way ice computation
41	REAL, DIMENSION(iip1,jjp1), INTENT(IN) :: masque ! land mask
42	#ifndef CPP_EARTH
43	WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
44	#else
45	!-------------------------------------------------------------------------------
46	! Local variables:
47	#include "control.h"
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	STOP
270
271
272	!===============================================================================
273	!
274	CONTAINS
275	!
276	!===============================================================================
277
278
279	!-------------------------------------------------------------------------------
280	!
281	SUBROUTINE get_2Dfield(fnam, mode, ibar, ndays, champo, flag, mask, lCPL)
282	!
283	!-------------------------------------------------------------------------------
284	! Comments:
285	! There are two assumptions concerning the NetCDF files, that are satisfied
286	! with files that are conforming NC convention:
287	! 1) The last dimension of the variables used is the time record.
288	! 2) Dimensional variables have the same names as corresponding dimensions.
289	!-------------------------------------------------------------------------------
290	USE netcdf, ONLY : NF90_OPEN, NF90_INQ_VARID, NF90_INQUIRE_VARIABLE, &
291	NF90_CLOSE, NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION, &
292	NF90_GET_VAR, NF90_GET_ATT
293	USE dimphy, ONLY : klon
294	USE phys_state_var_mod, ONLY : pctsrf
295	IMPLICIT NONE
296	#include "dimensions.h"
297	#include "paramet.h"
298	#include "comgeom2.h"
299	#include "control.h"
300	#include "indicesol.h"
301	#include "iniprint.h"
302	!-------------------------------------------------------------------------------
303	! Arguments:
304	CHARACTER(LEN=*), INTENT(IN) :: fnam ! NetCDF file name
305	CHARACTER(LEN=3), INTENT(IN) :: mode ! RUG, SIC, SST or ALB
306	LOGICAL, INTENT(IN) :: ibar ! interp on pressure levels
307	INTEGER, INTENT(IN) :: ndays ! current year number of days
308	REAL, POINTER, DIMENSION(:,:) :: champo ! output field = f(t)
309	LOGICAL, OPTIONAL, INTENT(IN) :: flag ! extrapol. (SST) old ice (SIC)
310	REAL, OPTIONAL, DIMENSION(iim,jjp1), INTENT(IN) :: mask
311	LOGICAL, OPTIONAL, INTENT(IN) :: lCPL ! Coupled model flag (for ICE)
312	!-------------------------------------------------------------------------------
313	! Local variables:
314	!--- NetCDF
315	INTEGER :: ierr, ncid, varid ! NetCDF identifiers
316	CHARACTER(LEN=30) :: dnam ! dimension name
317	CHARACTER(LEN=80) :: varname ! NetCDF variable name
318	!--- dimensions
319	INTEGER, DIMENSION(4) :: dids ! NetCDF dimensions identifiers
320	REAL, ALLOCATABLE, DIMENSION(:) :: dlon_ini ! initial longitudes vector
321	REAL, ALLOCATABLE, DIMENSION(:) :: dlat_ini ! initial latitudes vector
322	REAL, POINTER, DIMENSION(:) :: dlon, dlat ! reordered lon/lat vectors
323	!--- fields
324	INTEGER :: imdep, jmdep, lmdep ! dimensions of 'champ'
325	REAL, ALLOCATABLE, DIMENSION(:,:) :: champ ! wanted field on initial grid
326	REAL, ALLOCATABLE, DIMENSION(:) :: yder, timeyear
327	REAL, DIMENSION(iim,jjp1) :: champint ! interpolated field
328	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: champtime
329	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: champan
330	REAL :: time, by
331	!--- input files
332	CHARACTER(LEN=20) :: cal_in ! calendar
333	INTEGER :: ndays_in ! number of days
334	!--- misc
335	INTEGER :: i, j, k, l ! loop counters
336	REAL, ALLOCATABLE, DIMENSION(:,:) :: work ! used for extrapolation
337	CHARACTER(LEN=25) :: title ! for messages
338	LOGICAL :: extrp ! flag for extrapolation
339	REAL :: chmin, chmax
340	!-------------------------------------------------------------------------------
341	!---Variables depending on keyword 'mode' --------------------------------------
342	NULLIFY(champo)
343	SELECT CASE(mode)
344	CASE('RUG'); varname='RUGOS'; title='Rugosite'
345	CASE('SIC'); varname='sicbcs'; title='Sea-ice'
346	CASE('SST'); varname='tosbcs'; title='SST'
347	CASE('ALB'); varname='ALBEDO'; title='Albedo'
348	END SELECT
349	extrp=(flag.AND.mode=='SST')
350
351	!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE ------------------------------
352	ierr=NF90_OPEN(fnam,NF90_NOWRITE,ncid); CALL ncerr(ierr,fnam)
353	ierr=NF90_INQ_VARID(ncid,varname,varid); CALL ncerr(ierr,fnam)
354	ierr=NF90_INQUIRE_VARIABLE(ncid,varid,dimids=dids); CALL ncerr(ierr,fnam)
355
356	!--- Longitude
357	ierr=NF90_INQUIRE_DIMENSION(ncid,dids(1),name=dnam,len=imdep)
358	CALL ncerr(ierr,fnam); ALLOCATE(dlon_ini(imdep),dlon(imdep))
359	ierr=NF90_INQ_VARID(ncid,dnam,varid); CALL ncerr(ierr,fnam)
360	ierr=NF90_GET_VAR(ncid,varid,dlon_ini); CALL ncerr(ierr,fnam)
361	WRITE(lunout,*) 'variable ', dnam, 'dimension ', imdep
362
363	!--- Latitude
364	ierr=NF90_INQUIRE_DIMENSION(ncid,dids(2),name=dnam,len=jmdep)
365	CALL ncerr(ierr,fnam); ALLOCATE(dlat_ini(jmdep),dlat(jmdep))
366	ierr=NF90_INQ_VARID(ncid,dnam,varid); CALL ncerr(ierr,fnam)
367	ierr=NF90_GET_VAR(ncid,varid,dlat_ini); CALL ncerr(ierr,fnam)
368	WRITE(lunout,*) 'variable ', dnam, 'dimension ', jmdep
369
370	!--- Time (variable is not needed - it is rebuilt - but calendar is)
371	ierr=NF90_INQUIRE_DIMENSION(ncid,dids(3),name=dnam,len=lmdep)
372	CALL ncerr(ierr,fnam); ALLOCATE(timeyear(lmdep))
373	ierr=NF90_INQ_VARID(ncid,dnam,varid); CALL ncerr(ierr,fnam)
374	cal_in=' '
375	ierr=NF90_GET_ATT(ncid,varid,'calendar',cal_in)
376	IF(ierr/=NF90_NOERR) THEN
377	SELECT CASE(mode)
378	CASE('RUG','ALB'); cal_in='360d'
379	CASE('SIC','SST'); cal_in='gregorian'
380	END SELECT
381	WRITE(lunout,*)'ATTENTION: variable ''time'' sans attribut ''calendrier'' d&
382	&ans '//TRIM(fnam)//'. On choisit la valeur par defaut.'
383	END IF
384	WRITE(lunout,*) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', cal_in
385
386	!--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION ----------------------
387	!--- Determining input file number of days, depending on calendar
388	ndays_in=year_len(anneeref,cal_in)
389
390	!--- Time vector reconstruction (time vector from file is not trusted)
391	!--- If input records are not monthly, time sampling has to be constant !
392	timeyear=mid_months(anneeref,cal_in,lmdep)
393	IF(lmdep/=12) WRITE(lunout,'(a,i3,a)')'Note: les fichiers de '//TRIM(mode) &
394	//' ne comportent pas 12, mais ',lmdep,' enregistrements.'
395
396	!--- GETTING THE FIELD AND INTERPOLATING IT ------------------------------------
397	ALLOCATE(champ(imdep,jmdep),champtime(iim,jjp1,lmdep))
398	IF(extrp) ALLOCATE(work(imdep,jmdep))
399
400	WRITE(lunout,*)
401	WRITE(lunout,'(a,i3,a)')'LECTURE ET INTERPOLATION HORIZ. DE ',lmdep,' CHAMPS.'
402	ierr=NF90_INQ_VARID(ncid,varname,varid); CALL ncerr(ierr,fnam)
403	DO l=1,lmdep
404	ierr=NF90_GET_VAR(ncid,varid,champ,(/1,1,l/),(/imdep,jmdep,1/))
405	CALL ncerr(ierr,fnam)
406	CALL conf_dat2d(title,imdep,jmdep,dlon_ini,dlat_ini,dlon,dlat,champ,ibar)
407
408	IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
409
410	IF(ibar.AND..NOT.(mode=='SIC'.AND.flag)) THEN
411	IF(l==1) THEN
412	WRITE(lunout,*) &
413	'-------------------------------------------------------------------------'
414	WRITE(lunout,*) &
415	'$$$ Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$'
416	WRITE(lunout,*) &
417	'-------------------------------------------------------------------------'
418	END IF
419	IF(mode=='RUG') champ=LOG(champ)
420	CALL inter_barxy(dlon, dlat(:jmdep-1), champ, rlonu(:iim), &
421	rlatv, champint)
422	IF(mode=='RUG') THEN
423	champint=EXP(champint)
424	WHERE(NINT(mask)/=1) champint=0.001
425	END IF
426	ELSE
427	SELECT CASE(mode)
428	CASE('RUG'); CALL rugosite(imdep, jmdep, dlon, dlat, champ, &
429	iim, jjp1, rlonv, rlatu, champint, mask)
430	CASE('SIC'); CALL sea_ice (imdep, jmdep, dlon, dlat, champ, &
431	iim, jjp1, rlonv, rlatu, champint)
432	CASE('SST','ALB'); CALL grille_m(imdep, jmdep, dlon, dlat, champ, &
433	iim, jjp1, rlonv, rlatu, champint)
434	END SELECT
435	END IF
436	champtime(:,:,l)=champint
437	END DO
438	ierr=NF90_CLOSE(ncid); CALL ncerr(ierr,fnam)
439
440	DEALLOCATE(dlon_ini,dlat_ini,dlon,dlat,champ)
441	IF(extrp) DEALLOCATE(work)
442
443	!--- TIME INTERPOLATION --------------------------------------------------------
444	WRITE(lunout,*)
445	WRITE(lunout,*)'INTERPOLATION TEMPORELLE.'
446	WRITE(lunout,"(2x,' Vecteur temps en entree: ',10f6.1)") timeyear
447	WRITE(lunout,"(2x,' Vecteur temps en sortie de 0 a ',i3)") ndays
448	ALLOCATE(yder(lmdep),champan(iip1,jjp1,ndays))
449	DO j=1,jjp1
450	DO i=1,iim
451	CALL spline(timeyear,champtime(i,j,:),lmdep,1.e30,1.e30,yder)
452	DO k=1,ndays
453	time=FLOAT((k-1)*ndays_in)/FLOAT(ndays)
454	CALL splint(timeyear,champtime(i,j,:),yder,lmdep,time,by)
455	champan(i,j,k) = by
456	END DO
457	END DO
458	END DO
459	champan(iip1,:,:)=champan(1,:,:)
460	DEALLOCATE(yder,champtime,timeyear)
461
462	!--- Checking the result
463	DO j=1,jjp1
464	CALL minmax(iip1,champan(1,j,10),chmin,chmax)
465	WRITE(lunout,*)' '//TRIM(title)//' au temps 10 ',chmin,chmax,j
466	END DO
467
468	!--- SPECIAL FILTER FOR SST: SST>271.38 ----------------------------------------
469	IF(mode=='SST') THEN
470	WRITE(lunout,*) 'Filtrage de la SST: SST >= 271.38'
471	WHERE(champan<271.38) champan=271.38
472	END IF
473
474	!--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 ---------------------------------------
475	IF(mode=='SIC') THEN
476	WRITE(lunout,*) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
477	IF(.NOT.lCPL) champan(:,:,:)=champan(:,:,:)/100.
478	champan(iip1,:,:)=champan(1,:,:)
479	WHERE(champan>1.0) champan=1.0
480	WHERE(champan<0.0) champan=0.0
481	END IF
482
483	!--- DYNAMICAL TO PHYSICAL GRID ------------------------------------------------
484	ALLOCATE(champo(klon,ndays))
485	DO k=1,ndays
486	CALL gr_dyn_fi(1,iip1,jjp1,klon,champan(1,1,k),champo(1,k))
487	END DO
488	DEALLOCATE(champan)
489
490	END SUBROUTINE get_2Dfield
491	!
492	!-------------------------------------------------------------------------------
493
494
495
496	!-------------------------------------------------------------------------------
497	!
498	FUNCTION year_len(y,cal_in)
499	!
500	!-------------------------------------------------------------------------------
501	USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_year_len
502	IMPLICIT NONE
503	!-------------------------------------------------------------------------------
504	! Arguments:
505	INTEGER :: year_len
506	INTEGER, INTENT(IN) :: y
507	CHARACTER(LEN=*), INTENT(IN) :: cal_in
508	!-------------------------------------------------------------------------------
509	! Local variables:
510	CHARACTER(LEN=20) :: cal_out ! calendar (for outputs)
511	!-------------------------------------------------------------------------------
512	!--- Getting the input calendar to reset at the end of the function
513	CALL ioget_calendar(cal_out)
514
515	!--- Unlocking calendar and setting it to wanted one
516	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
517
518	!--- Getting the number of days in this year
519	year_len=ioget_year_len(y)
520
521	!--- Back to original calendar
522	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
523
524	END FUNCTION year_len
525	!
526	!-------------------------------------------------------------------------------
527
528
529	!-------------------------------------------------------------------------------
530	!
531	FUNCTION mid_months(y,cal_in,nm)
532	!
533	!-------------------------------------------------------------------------------
534	USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_mon_len
535	IMPLICIT NONE
536	!-------------------------------------------------------------------------------
537	! Arguments:
538	INTEGER, INTENT(IN) :: y ! year
539	CHARACTER(LEN=*), INTENT(IN) :: cal_in ! calendar
540	INTEGER, INTENT(IN) :: nm ! months/year number
541	REAL, DIMENSION(nm) :: mid_months ! mid-month times
542	!-------------------------------------------------------------------------------
543	! Local variables:
544	CHARACTER(LEN=99) :: mess ! error message
545	CHARACTER(LEN=20) :: cal_out ! calendar (for outputs)
546	INTEGER, DIMENSION(nm) :: mnth ! months lengths (days)
547	INTEGER :: m ! months counter
548	INTEGER :: nd ! number of days
549	!-------------------------------------------------------------------------------
550	nd=year_len(y,cal_in)
551
552	IF(nm==12) THEN
553
554	!--- Getting the input calendar to reset at the end of the function
555	CALL ioget_calendar(cal_out)
556
557	!--- Unlocking calendar and setting it to wanted one
558	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
559
560	!--- Getting the length of each month
561	DO m=1,nm; mnth(m)=ioget_mon_len(y,m); END DO
562
563	!--- Back to original calendar
564	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
565
566	ELSE IF(MODULO(nd,nm)/=0) THEN
567	WRITE(mess,'(a,i3,a,i3,a)')'Unconsistent calendar: ',nd,' days/year, but ',&
568	nm,' months/year. Months number should divide days number.'
569	CALL abort_gcm('mid_months',TRIM(mess),1)
570
571	ELSE
572	mnth=(/(m,m=1,nm,nd/nm)/)
573	END IF
574
575	!--- Mid-months times
576	mid_months(1)=0.5*FLOAT(mnth(1))
577	DO k=2,nm
578	mid_months(k)=mid_months(k-1)+0.5*FLOAT(mnth(k-1)+mnth(k))
579	END DO
580
581	END FUNCTION mid_months
582	!
583	!-------------------------------------------------------------------------------
584
585
586	!-------------------------------------------------------------------------------
587	!
588	SUBROUTINE ncerr(ncres,fnam)
589	!
590	!-------------------------------------------------------------------------------
591	! Purpose: NetCDF errors handling.
592	!-------------------------------------------------------------------------------
593	USE netcdf, ONLY : NF90_NOERR, NF90_STRERROR
594	IMPLICIT NONE
595	!-------------------------------------------------------------------------------
596	! Arguments:
597	INTEGER, INTENT(IN) :: ncres
598	CHARACTER(LEN=*), INTENT(IN) :: fnam
599	!-------------------------------------------------------------------------------
600	#include "iniprint.h"
601	IF(ncres/=NF90_NOERR) THEN
602	WRITE(lunout,*)'Problem with file '//TRIM(fnam)//' in routine limit_netcdf.'
603	CALL abort_gcm('limit_netcdf',NF90_STRERROR(ncres),1)
604	END IF
605
606	END SUBROUTINE ncerr
607	!
608	!-------------------------------------------------------------------------------
609
610
611	END SUBROUTINE limit_netcdf

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

Download in other formats:

Original Format