Context Navigation

limit_netcdf.F90 @ 1969

Last change on this file since 1969 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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: 26.4 KB

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

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

Download in other formats:

Original Format