Context Navigation

limit_netcdf.F90 @ 2328

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

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

Download in other formats:

Original Format