Context Navigation

limit_netcdf.F90 @ 2310

Last change on this file since 2310 was 2302, checked in by Ehouarn Millour, 9 years ago
Move etat0phys_netcdf.F90 to "dynlonlat_phylonlat/phylmd" as it relies on "phylmd" routines. Some cleanup to remove obsolete and unecessary CPP_EARTH preprocessing condition. EM
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.2 KB

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

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

Download in other formats:

Original Format