Context Navigation

limit_netcdf.F90 @ 2311

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

Further modifications to enforce physics/dynamics separation:

moved iniprint.h and misc_mod back to dyn3d_common, as these should only be used by dynamics.
created print_control_mod in the physics to store flags prt_level, lunout, debug to be local to physics (should be used rather than iniprint.h)
created abort_physic.F90 , which does the same job as abort_gcm() did, but should be used instead when in physics.
reactivated inifis (turned it into a module, inifis_mod.F90) to initialize physical constants and print_control_mod flags.

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.3 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 "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_physic('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_physic('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	!-----------------------------------------------------------------------------
280	! Arguments:
281	CHARACTER(LEN=*), INTENT(IN) :: fnam ! NetCDF file name
282	CHARACTER(LEN=10), INTENT(IN) :: varname ! NetCDF variable name
283	CHARACTER(LEN=3), INTENT(IN) :: mode ! RUG, SIC, SST or ALB
284	LOGICAL, INTENT(IN) :: ibar ! interp on pressure levels
285	INTEGER, INTENT(IN) :: ndays ! current year number of days
286	REAL, POINTER, DIMENSION(:, :) :: champo ! output field = f(t)
287	LOGICAL, OPTIONAL, INTENT(IN) :: flag ! extrapol. (SST) old ice (SIC)
288	REAL, OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask
289	!------------------------------------------------------------------------------
290	! Local variables:
291	!--- NetCDF
292	INTEGER :: ncid, varid ! NetCDF identifiers
293	CHARACTER(LEN=30) :: dnam ! dimension name
294	!--- dimensions
295	INTEGER :: dids(4) ! NetCDF dimensions identifiers
296	REAL, ALLOCATABLE :: dlon_ini(:) ! initial longitudes vector
297	REAL, ALLOCATABLE :: dlat_ini(:) ! initial latitudes vector
298	REAL, POINTER :: dlon(:), dlat(:) ! reordered lon/lat vectors
299	!--- fields
300	INTEGER :: imdep, jmdep, lmdep ! dimensions of 'champ'
301	REAL, ALLOCATABLE :: champ(:,:) ! wanted field on initial grid
302	REAL, ALLOCATABLE :: yder(:), timeyear(:)
303	REAL :: champint(iim,jjp1) ! interpolated field
304	REAL, ALLOCATABLE :: champtime(:,:,:)
305	REAL, ALLOCATABLE :: champan(:,:,:)
306	!--- input files
307	CHARACTER(LEN=20) :: cal_in ! calendar
308	CHARACTER(LEN=20) :: unit_sic ! attribute unit in sea-ice file
309	INTEGER :: ndays_in ! number of days
310	!--- misc
311	INTEGER :: i, j, k, l ! loop counters
312	REAL, ALLOCATABLE :: work(:,:) ! used for extrapolation
313	CHARACTER(LEN=25) :: title ! for messages
314	LOGICAL :: extrp ! flag for extrapolation
315	LOGICAL :: oldice ! flag for old way ice computation
316	REAL :: chmin, chmax
317	INTEGER ierr
318	integer n_extrap ! number of extrapolated points
319	logical skip
320
321	!------------------------------------------------------------------------------
322	!---Variables depending on keyword 'mode' -------------------------------------
323	NULLIFY(champo)
324
325	SELECT CASE(mode)
326	CASE('RUG'); title='Rugosite'
327	CASE('SIC'); title='Sea-ice'
328	CASE('SST'); title='SST'
329	CASE('ALB'); title='Albedo'
330	END SELECT
331
332
333	extrp=.FALSE.
334	oldice=.FALSE.
335	IF ( PRESENT(flag) ) THEN
336	IF ( flag .AND. mode=='SST' ) extrp=.TRUE.
337	IF ( flag .AND. mode=='SIC' ) oldice=.TRUE.
338	END IF
339
340	!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
341	CALL msg(5,' Now reading file : '//TRIM(fnam))
342	CALL ncerr(NF90_OPEN(fnam, NF90_NOWRITE, ncid),fnam)
343	CALL ncerr(NF90_INQ_VARID(ncid, trim(varname), varid),fnam)
344	CALL ncerr(NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids),fnam)
345
346	!--- Read unit for sea-ice variable only
347	IF (mode=='SIC') THEN
348	IF(NF90_GET_ATT(ncid, varid, 'units', unit_sic)/=NF90_NOERR) THEN
349	CALL msg(5,'No unit in sea-ice file. Take percentage as default value')
350	unit_sic='X'
351	ELSE
352	CALL msg(5,'Sea-ice cover has unit='//TRIM(unit_sic))
353	END IF
354	END IF
355
356	!--- Longitude
357	CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(1), name=dnam, len=imdep),fnam)
358	ALLOCATE(dlon_ini(imdep), dlon(imdep))
359	CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
360	CALL ncerr(NF90_GET_VAR(ncid, varid, dlon_ini), fnam)
361	CALL msg(5,'variable '//TRIM(dnam)//' dimension ', imdep)
362
363	!--- Latitude
364	CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(2), name=dnam, len=jmdep),fnam)
365	ALLOCATE(dlat_ini(jmdep), dlat(jmdep))
366	CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
367	CALL ncerr(NF90_GET_VAR(ncid, varid, dlat_ini), fnam)
368	CALL msg(5,'variable '//TRIM(dnam)//' dimension ', jmdep)
369
370	!--- Time (variable is not needed - it is rebuilt - but calendar is)
371	CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(3), name=dnam, len=lmdep), fnam)
372	ALLOCATE(timeyear(lmdep))
373	CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
374	cal_in=' '
375	IF(NF90_GET_ATT(ncid, varid, 'calendar', cal_in)/=NF90_NOERR) THEN
376	SELECT CASE(mode)
377	CASE('RUG', 'ALB'); cal_in='360d'
378	CASE('SIC', 'SST'); cal_in='gregorian'
379	END SELECT
380	CALL msg(5,'WARNING: missing "calendar" attribute for "time" in '&
381	&//TRIM(fnam)//'. Choosing default value.')
382	END IF
383	CALL msg(5,'var, calendar, dim: '//TRIM(dnam)//' '//TRIM(cal_in), lmdep)
384
385	!--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION --------------------
386	!--- Determining input file number of days, depending on calendar
387	ndays_in=year_len(anneeref, cal_in)
388
389	!--- Time vector reconstruction (time vector from file is not trusted)
390	!--- If input records are not monthly, time sampling has to be constant !
391	timeyear=mid_months(anneeref, cal_in, lmdep)
392	IF (lmdep /= 12) WRITE(lunout,*) 'Note : les fichiers de ', TRIM(mode), &
393	' ne comportent pas 12, mais ', lmdep, ' enregistrements.'
394
395	!--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
396	ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep))
397	IF(extrp) ALLOCATE(work(imdep, jmdep))
398	CALL msg(5,'')
399	CALL msg(5,'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, ' CHAMPS.')
400	CALL ncerr(NF90_INQ_VARID(ncid, varname, varid), fnam)
401	DO l=1, lmdep
402	CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,l],[imdep,jmdep,1]),fnam)
403	CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, ibar)
404	IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
405
406	IF(ibar .AND. .NOT.oldice) THEN
407	IF(l==1) THEN
408	CALL msg(5,"----------------------------------------------------------")
409	CALL msg(5,"$$$ Interpolation barycentrique pour "//TRIM(title)//" $$$")
410	CALL msg(5,"----------------------------------------------------------")
411	END IF
412	IF(mode=='RUG') champ=LOG(champ)
413	CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint)
414	IF(mode=='RUG') THEN
415	champint=EXP(champint)
416	WHERE(NINT(mask)/=1) champint=0.001
417	END IF
418	ELSE
419	SELECT CASE(mode)
420	CASE('RUG'); CALL rugosite(dlon,dlat,champ,rlonv,rlatu,champint,mask)
421	CASE('SIC'); CALL sea_ice (dlon,dlat,champ,rlonv,rlatu,champint)
422	CASE('SST','ALB'); CALL grille_m(dlon,dlat,champ,rlonv,rlatu,champint)
423	END SELECT
424	END IF
425	champtime(:, :, l)=champint
426	END DO
427	CALL ncerr(NF90_CLOSE(ncid), fnam)
428
429	DEALLOCATE(dlon_ini, dlat_ini, dlon, dlat, champ)
430	IF(extrp) DEALLOCATE(work)
431
432	!--- TIME INTERPOLATION ------------------------------------------------------
433	IF(prt_level>5) THEN
434	WRITE(lunout, *)
435	WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
436	WRITE(lunout, *)' Vecteur temps en entree: ', timeyear
437	WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays
438	END IF
439
440	ALLOCATE(yder(lmdep), champan(iip1, jjp1, ndays))
441	skip = .false.
442	n_extrap = 0
443	DO j=1, jjp1
444	DO i=1, iim
445	yder = pchsp_95(timeyear, champtime(i, j, :), ibeg=2, iend=2, &
446	vc_beg=0., vc_end=0.)
447	CALL pchfe_95(timeyear, champtime(i, j, :), yder, skip, &
448	arth(0., real(ndays_in) / ndays, ndays), champan(i, j, :), ierr)
449	if (ierr < 0) stop 1
450	n_extrap = n_extrap + ierr
451	END DO
452	END DO
453	if (n_extrap /= 0) then
454	WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
455	end if
456	champan(iip1, :, :)=champan(1, :, :)
457	DEALLOCATE(yder, champtime, timeyear)
458
459	!--- Checking the result
460	DO j=1, jjp1
461	CALL minmax(iip1, champan(1, j, 10), chmin, chmax)
462	IF (prt_level>5) WRITE(lunout, *)' ',TRIM(title),' au temps 10 ', chmin, chmax, j
463	END DO
464
465	!--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
466	IF(mode=='SST') THEN
467	CALL msg(5,'Filtrage de la SST: SST >= 271.38')
468	WHERE(champan<271.38) champan=271.38
469	END IF
470
471	!--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
472	IF(mode=='SIC') THEN
473	CALL msg(5,'Filtrage de la SIC: 0.0 < Sea-ice < 1.0')
474
475	IF (unit_sic=='1') THEN
476	! Nothing to be done for sea-ice field is already in fraction of 1
477	! This is the case for sea-ice in file cpl_atm_sic.nc
478	CALL msg(5,'Sea-ice field already in fraction of 1')
479	ELSE
480	! Convert sea ice from percentage to fraction of 1
481	CALL msg(5,'Transformt sea-ice field from percentage to fraction of 1.')
482	champan(:, :, :)=champan(:, :, :)/100.
483	END IF
484
485	champan(iip1, :, :)=champan(1, :, :)
486	WHERE(champan>1.0) champan=1.0
487	WHERE(champan<0.0) champan=0.0
488	END IF
489
490	!--- DYNAMICAL TO PHYSICAL GRID ----------------------------------------------
491	ALLOCATE(champo(klon, ndays))
492	DO k=1, ndays
493	CALL gr_dyn_fi(1, iip1, jjp1, klon, champan(1, 1, k), champo(1, k))
494	END DO
495	DEALLOCATE(champan)
496
497	END SUBROUTINE get_2Dfield
498	!
499	!-------------------------------------------------------------------------------
500
501
502	!-------------------------------------------------------------------------------
503	!
504	FUNCTION year_len(y,cal_in)
505	!
506	!-------------------------------------------------------------------------------
507	USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_year_len
508	IMPLICIT NONE
509	!-------------------------------------------------------------------------------
510	! Arguments:
511	INTEGER :: year_len
512	INTEGER, INTENT(IN) :: y
513	CHARACTER(LEN=*), INTENT(IN) :: cal_in
514	!-------------------------------------------------------------------------------
515	! Local variables:
516	CHARACTER(LEN=20) :: cal_out ! calendar (for outputs)
517	!-------------------------------------------------------------------------------
518	!--- Getting the input calendar to reset at the end of the function
519	CALL ioget_calendar(cal_out)
520
521	!--- Unlocking calendar and setting it to wanted one
522	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
523
524	!--- Getting the number of days in this year
525	year_len=ioget_year_len(y)
526
527	!--- Back to original calendar
528	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
529
530	END FUNCTION year_len
531	!
532	!-------------------------------------------------------------------------------
533
534
535	!-------------------------------------------------------------------------------
536	!
537	FUNCTION mid_months(y,cal_in,nm)
538	!
539	!-------------------------------------------------------------------------------
540	USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_mon_len
541	IMPLICIT NONE
542	!-------------------------------------------------------------------------------
543	! Arguments:
544	INTEGER, INTENT(IN) :: y ! year
545	CHARACTER(LEN=*), INTENT(IN) :: cal_in ! calendar
546	INTEGER, INTENT(IN) :: nm ! months/year number
547	REAL, DIMENSION(nm) :: mid_months ! mid-month times
548	!-------------------------------------------------------------------------------
549	! Local variables:
550	CHARACTER(LEN=99) :: mess ! error message
551	CHARACTER(LEN=20) :: cal_out ! calendar (for outputs)
552	INTEGER, DIMENSION(nm) :: mnth ! months lengths (days)
553	INTEGER :: m ! months counter
554	INTEGER :: nd ! number of days
555	!-------------------------------------------------------------------------------
556	nd=year_len(y,cal_in)
557
558	IF(nm==12) THEN
559
560	!--- Getting the input calendar to reset at the end of the function
561	CALL ioget_calendar(cal_out)
562
563	!--- Unlocking calendar and setting it to wanted one
564	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
565
566	!--- Getting the length of each month
567	DO m=1,nm; mnth(m)=ioget_mon_len(y,m); END DO
568
569	!--- Back to original calendar
570	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
571
572	ELSE IF(MODULO(nd,nm)/=0) THEN
573	WRITE(mess,'(a,i3,a,i3,a)')'Unconsistent calendar: ',nd,' days/year, but ',&
574	nm,' months/year. Months number should divide days number.'
575	CALL abort_physic('mid_months',TRIM(mess),1)
576
577	ELSE
578	mnth=[(m,m=1,nm,nd/nm)]
579	END IF
580
581	!--- Mid-months times
582	mid_months(1)=0.5*REAL(mnth(1))
583	DO k=2,nm
584	mid_months(k)=mid_months(k-1)+0.5*REAL(mnth(k-1)+mnth(k))
585	END DO
586
587	END FUNCTION mid_months
588	!
589	!-------------------------------------------------------------------------------
590
591
592
593	!-------------------------------------------------------------------------------
594	!
595	SUBROUTINE msg(lev,str1,i,str2)
596	!
597	!-------------------------------------------------------------------------------
598	! Arguments:
599	INTEGER, INTENT(IN) :: lev
600	CHARACTER(LEN=*), INTENT(IN) :: str1
601	INTEGER, OPTIONAL, INTENT(IN) :: i
602	CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: str2
603	!-------------------------------------------------------------------------------
604	IF(prt_level>lev) THEN
605	IF(PRESENT(str2)) THEN
606	WRITE(lunout,*) TRIM(str1), i, TRIM(str2)
607	ELSE IF(PRESENT(i)) THEN
608	WRITE(lunout,*) TRIM(str1), i
609	ELSE
610	WRITE(lunout,*) TRIM(str1)
611	END IF
612	END IF
613
614	END SUBROUTINE msg
615	!
616	!-------------------------------------------------------------------------------
617
618
619	!-------------------------------------------------------------------------------
620	!
621	SUBROUTINE ncerr(ncres,fnam)
622	!
623	!-------------------------------------------------------------------------------
624	! Purpose: NetCDF errors handling.
625	!-------------------------------------------------------------------------------
626	USE netcdf, ONLY : NF90_NOERR, NF90_STRERROR
627	USE print_control_mod, ONLY: lunout
628	IMPLICIT NONE
629	!-------------------------------------------------------------------------------
630	! Arguments:
631	INTEGER, INTENT(IN) :: ncres
632	CHARACTER(LEN=*), INTENT(IN) :: fnam
633	!-------------------------------------------------------------------------------
634	IF(ncres/=NF90_NOERR) THEN
635	WRITE(lunout,*)'Problem with file '//TRIM(fnam)//' in routine limit_netcdf.'
636	CALL abort_physic('limit_netcdf',NF90_STRERROR(ncres),1)
637	END IF
638
639	END SUBROUTINE ncerr
640	!
641	!-------------------------------------------------------------------------------
642
643	#endif
644	! of #ifndef CPP_1D
645	END SUBROUTINE limit_netcdf

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

Download in other formats:

Original Format