Context Navigation

limit_netcdf.F90 @ 2239

Last change on this file since 2239 was 2239, checked in by Ehouarn Millour, 9 years ago

Reorganizing physics/dynamics interface:

what is related to dynamics-physics interface is now in a seperate directory: dynlmdz_phy* for physics in phy*
1d model and related dependencies (including a couple from "dynamics", set up as symbolic links) is now in subdirectory "dyn1d" of phy*.
"bibio" directory is now "misc" and should only contain autonomous utilities.
"cosp" is now a subdirectory of phylmd.

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

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

Download in other formats:

Original Format