Context Navigation

limit_netcdf.F90 @ 2061

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

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

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

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

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

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

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

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

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

File size: 26.0 KB

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

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

Download in other formats:

Original Format