Context Navigation

limit_netcdf.F90 @ 2588

Last change on this file since 2588 was 2576, checked in by dcugnet, 8 years ago

Use "traditional" sst/sic files without additional records from contiguous years,
but check for previous and next files. Example for Amip Sea Surface Temperature:

amipbc_sst_1x1.nc must be present.
if amipbc_sst_1x1_m.nc is present (previous year), its last record is used

for interpolation ; otherwise, the last record of amipbc_sst_1x1.nc is used.

if amipbc_sst_1x1_p.nc is present (next year), its first record is used

for interpolation ; otherwise, the first record of amipbc_sst_1x1.nc is used.
If both contiguous fiels are missing, the interpolation then implies periodicity.
Note: grid_noro0 has been moved from limit_netcd to grid_noro_m where it belongs.

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

Line
1	MODULE limit
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
19	USE ioipsl, ONLY: flininfo, flinopen, flinget, flinclo, &
20	ioconf_calendar, ioget_calendar, lock_calendar, ioget_mon_len, ioget_year_len
21	USE assert_eq_m, ONLY: assert_eq
22	USE conf_dat_m, ONLY: conf_dat2d, conf_dat3d
23	USE dimphy, ONLY: klon, zmasq
24	USE geometry_mod, ONLY: longitude_deg, latitude_deg
25	USE phys_state_var_mod, ONLY: pctsrf
26	USE control_mod, ONLY: anneeref
27	USE init_ssrf_m, ONLY: start_init_subsurf
28
29	CHARACTER(LEN=20), PARAMETER :: &
30	fsst(4)=['amipbc_sst_1x1.nc ','cpl_atm_sst.nc ','histmth_sst.nc '&
31	,'sstk.nc ']
32	CHARACTER(LEN=20), PARAMETER :: &
33	fsic(4)=['amipbc_sic_1x1.nc ','cpl_atm_sic.nc ','histmth_sic.nc '&
34	,'ci.nc ']
35	CHARACTER(LEN=10), PARAMETER :: &
36	vsst(4)=['tosbcs ','SISUTESW ','tsol_oce ','sstk '], &
37	vsic(4)=['sicbcs ','SIICECOV ','pourc_sic ','ci ']
38	CHARACTER(LEN=10), PARAMETER :: &
39	frugo='Rugos.nc ', falbe='Albedo.nc ', frelf='Relief.nc ', &
40	vrug='RUGOS ', valb='ALBEDO ', vrel='RELIEF '
41
42	CONTAINS
43
44	!-------------------------------------------------------------------------------
45	!
46	SUBROUTINE limit_netcdf(masque, phis, extrap)
47	!
48	!-------------------------------------------------------------------------------
49	! Author : L. Fairhead, 27/01/94
50	!-------------------------------------------------------------------------------
51	! Purpose: Boundary conditions files building for new model using climatologies.
52	! Both grids have to be regular.
53	!-------------------------------------------------------------------------------
54	! Note: This routine is designed to work for Earth
55	!-------------------------------------------------------------------------------
56	! Modification history:
57	! * 23/03/1994: Z. X. Li
58	! * 09/1999: L. Fairhead (netcdf reading in LMDZ.3.3)
59	! * 07/2001: P. Le Van
60	! * 11/2009: L. Guez (ozone day & night climatos, see etat0_netcdf.F90)
61	! * 12/2009: D. Cugnet (f77->f90, calendars, files from coupled runs)
62	! * 04/2016: D. Cugnet (12/14 recs SST/SIC files: cyclic/interannual runs)
63	!-------------------------------------------------------------------------------
64	#ifndef CPP_1D
65	USE indice_sol_mod
66	USE netcdf, ONLY: NF90_OPEN, NF90_CREATE, NF90_CLOSE, &
67	NF90_DEF_DIM, NF90_DEF_VAR, NF90_PUT_VAR, NF90_PUT_ATT, &
68	NF90_NOERR, NF90_NOWRITE, NF90_DOUBLE, NF90_GLOBAL, &
69	NF90_CLOBBER, NF90_ENDDEF, NF90_UNLIMITED, NF90_FLOAT
70	USE inter_barxy_m, ONLY: inter_barxy
71	USE netcdf95, ONLY: nf95_def_var, nf95_put_att, nf95_put_var
72	IMPLICIT NONE
73	!-------------------------------------------------------------------------------
74	! Arguments:
75	include "iniprint.h"
76	include "dimensions.h"
77	include "paramet.h"
78	REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: masque ! land mask
79	REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: phis ! ground geopotential
80	LOGICAL, INTENT(IN) :: extrap ! SST extrapolation flag
81	!-------------------------------------------------------------------------------
82	! Local variables:
83	include "logic.h"
84	include "comgeom2.h"
85	include "comconst.h"
86
87	!--- INPUT NETCDF FILES NAMES --------------------------------------------------
88	CHARACTER(LEN=20) :: icefile, sstfile, dumstr, fnam
89	CHARACTER(LEN=10) :: varname
90
91	!--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------
92	REAL :: fi_ice(klon), verif(klon)
93	REAL, POINTER :: phy_rug(:,:)=>NULL(), phy_ice(:,:)=>NULL()
94	REAL, POINTER :: phy_sst(:,:)=>NULL(), phy_alb(:,:)=>NULL()
95	REAL, ALLOCATABLE :: phy_bil(:,:), pctsrf_t(:,:,:)
96	INTEGER :: nbad
97
98	!--- VARIABLES FOR OUTPUT FILE WRITING -----------------------------------------
99	INTEGER :: ierr, nid, ndim, ntim, k, dims(2), ix_sic, ix_sst
100	INTEGER :: id_tim, id_SST, id_BILS, id_RUG, id_ALB
101	INTEGER :: id_FOCE, id_FSIC, id_FTER, id_FLIC, varid_longitude, varid_latitude
102	INTEGER :: NF90_FORMAT
103	INTEGER :: ndays !--- Depending on the output calendar
104
105	!--- INITIALIZATIONS -----------------------------------------------------------
106	#ifdef NC_DOUBLE
107	NF90_FORMAT=NF90_DOUBLE
108	#else
109	NF90_FORMAT=NF90_FLOAT
110	#endif
111	CALL inigeom
112
113	!--- MASK, GROUND GEOPOT. & SUBSURFACES COMPUTATION (IN CASE ok_etat0==.FALSE.)
114	IF(ALL(masque==-99999.)) THEN
115	CALL start_init_orog0(rlonv,rlatu,phis,masque)
116	CALL gr_dyn_fi(1,iip1,jjp1,klon,masque,zmasq) !--- To physical grid
117	ALLOCATE(pctsrf(klon,nbsrf))
118	CALL start_init_subsurf(.FALSE.)
119	!--- TO MATCH EXACTLY WHAT WOULD BE DONE IN etat0phys_netcdf
120	WHERE( masque(:,:)<EPSFRA) masque(:,:)=0.
121	WHERE(1.-masque(:,:)<EPSFRA) masque(:,:)=1.
122	END IF
123
124	!--- Beware: anneeref (from gcm.def) is used to determine output time sampling
125	ndays=ioget_year_len(anneeref)
126
127	!--- RUGOSITY TREATMENT --------------------------------------------------------
128	CALL msg(1,'Traitement de la rugosite')
129	CALL get_2Dfield(frugo,vrug,'RUG',ndays,phy_rug,mask=masque(1:iim,:))
130
131	!--- OCEAN TREATMENT -----------------------------------------------------------
132	CALL msg(1,'Traitement de la glace oceanique')
133
134	! Input SIC file selection
135	! Open file only to test if available
136	DO ix_sic=1,SIZE(fsic)
137	IF ( NF90_OPEN(TRIM(fsic(ix_sic)),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
138	icefile=fsic(ix_sic); varname=vsic(ix_sic); EXIT
139	END IF
140	END DO
141	IF(ix_sic==SIZE(fsic)+1) THEN
142	WRITE(lunout,*) 'ERROR! No sea-ice input file was found.'
143	WRITE(lunout,*) 'One of following files must be available : '
144	DO k=1,SIZE(fsic); WRITE(lunout,*) TRIM(fsic(k)); END DO
145	CALL abort_physic('limit_netcdf','No sea-ice file was found',1)
146	END IF
147	CALL ncerr(NF90_CLOSE(nid),icefile)
148	CALL msg(-1,'Fichier choisi pour la glace de mer:'//TRIM(icefile))
149
150	CALL get_2Dfield(icefile,varname, 'SIC',ndays,phy_ice)
151
152	ALLOCATE(pctsrf_t(klon,nbsrf,ndays))
153	DO k=1,ndays
154	fi_ice=phy_ice(:,k)
155	WHERE(fi_ice>=1.0 ) fi_ice=1.0
156	WHERE(fi_ice<EPSFRA) fi_ice=0.0
157	pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter) ! land soil
158	pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic) ! land ice
159	SELECT CASE(ix_sic)
160	CASE(2) ! SIC=pICE*(1-LIC-TER) (CPL)
161	pctsrf_t(:,is_sic,k)=fi_ice(:)*(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter))
162	CASE(3) ! SIC=pICE (HIST)
163	pctsrf_t(:,is_sic,k)=fi_ice(:)
164	CASE DEFAULT ! SIC=pICE-LIC (AMIP,ERAI)
165	pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k)
166	END SELECT
167	WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0.
168	WHERE(1.0-zmasq<EPSFRA)
169	pctsrf_t(:,is_sic,k)=0.0
170	pctsrf_t(:,is_oce,k)=0.0
171	ELSEWHERE
172	WHERE(pctsrf_t(:,is_sic,k)>=1.0-zmasq)
173	pctsrf_t(:,is_sic,k)=1.0-zmasq
174	pctsrf_t(:,is_oce,k)=0.0
175	ELSEWHERE
176	pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)
177	WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)
178	pctsrf_t(:,is_oce,k)=0.0
179	pctsrf_t(:,is_sic,k)=1.0-zmasq
180	END WHERE
181	END WHERE
182	END WHERE
183	nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0)
184	IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb points = ',nbad
185	nbad=COUNT(ABS(SUM(pctsrf_t(:,:,k),DIM=2)-1.0)>EPSFRA)
186	IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad
187	END DO
188	DEALLOCATE(phy_ice)
189
190	!--- SST TREATMENT -------------------------------------------------------------
191	CALL msg(1,'Traitement de la sst')
192
193	! Input SST file selection
194	! Open file only to test if available
195	DO ix_sst=1,SIZE(fsst)
196	IF ( NF90_OPEN(TRIM(fsst(ix_sst)),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
197	sstfile=fsst(ix_sst); varname=vsst(ix_sst); EXIT
198	END IF
199	END DO
200	IF(ix_sst==SIZE(fsst)+1) THEN
201	WRITE(lunout,*) 'ERROR! No sst input file was found.'
202	WRITE(lunout,*) 'One of following files must be available : '
203	DO k=1,SIZE(fsst); WRITE(lunout,*) TRIM(fsst(k)); END DO
204	CALL abort_physic('limit_netcdf','No sst file was found',1)
205	END IF
206	CALL ncerr(NF90_CLOSE(nid),sstfile)
207	CALL msg(-1,'Fichier choisi pour la temperature de mer: '//TRIM(sstfile))
208
209	CALL get_2Dfield(sstfile,varname,'SST',ndays,phy_sst,flag=extrap)
210
211	!--- ALBEDO TREATMENT ----------------------------------------------------------
212	CALL msg(1,'Traitement de l albedo')
213	CALL get_2Dfield(falbe,valb,'ALB',ndays,phy_alb)
214
215	!--- REFERENCE GROUND HEAT FLUX TREATMENT --------------------------------------
216	ALLOCATE(phy_bil(klon,ndays)); phy_bil=0.0
217
218	!--- OUTPUT FILE WRITING -------------------------------------------------------
219	CALL msg(5,'Ecriture du fichier limit : debut')
220	fnam="limit.nc"
221
222	!--- File creation
223	CALL ncerr(NF90_CREATE(fnam,NF90_CLOBBER,nid),fnam)
224	CALL ncerr(NF90_PUT_ATT(nid,NF90_GLOBAL,"title","Fichier conditions aux limites"),fnam)
225
226	!--- Dimensions creation
227	CALL ncerr(NF90_DEF_DIM(nid,"points_physiques",klon,ndim),fnam)
228	CALL ncerr(NF90_DEF_DIM(nid,"time",NF90_UNLIMITED,ntim),fnam)
229
230	dims=[ndim,ntim]
231
232	!--- Variables creation
233	CALL ncerr(NF90_DEF_VAR(nid,"TEMPS",NF90_FORMAT,[ntim],id_tim),fnam)
234	CALL ncerr(NF90_DEF_VAR(nid,"FOCE", NF90_FORMAT,dims,id_FOCE),fnam)
235	CALL ncerr(NF90_DEF_VAR(nid,"FSIC", NF90_FORMAT,dims,id_FSIC),fnam)
236	CALL ncerr(NF90_DEF_VAR(nid,"FTER", NF90_FORMAT,dims,id_FTER),fnam)
237	CALL ncerr(NF90_DEF_VAR(nid,"FLIC", NF90_FORMAT,dims,id_FLIC),fnam)
238	CALL ncerr(NF90_DEF_VAR(nid,"SST", NF90_FORMAT,dims,id_SST),fnam)
239	CALL ncerr(NF90_DEF_VAR(nid,"BILS", NF90_FORMAT,dims,id_BILS),fnam)
240	CALL ncerr(NF90_DEF_VAR(nid,"ALB", NF90_FORMAT,dims,id_ALB),fnam)
241	CALL ncerr(NF90_DEF_VAR(nid,"RUG", NF90_FORMAT,dims,id_RUG),fnam)
242	call nf95_def_var(nid, "longitude", NF90_FLOAT, ndim, varid_longitude)
243	call nf95_def_var(nid, "latitude", NF90_FLOAT, ndim, varid_latitude)
244
245	!--- Attributes creation
246	CALL ncerr(NF90_PUT_ATT(nid,id_tim, "title","Jour dans l annee"),fnam)
247	CALL ncerr(NF90_PUT_ATT(nid,id_FOCE,"title","Fraction ocean"),fnam)
248	CALL ncerr(NF90_PUT_ATT(nid,id_FSIC,"title","Fraction glace de mer"),fnam)
249	CALL ncerr(NF90_PUT_ATT(nid,id_FTER,"title","Fraction terre"),fnam)
250	CALL ncerr(NF90_PUT_ATT(nid,id_FLIC,"title","Fraction land ice"),fnam)
251	CALL ncerr(NF90_PUT_ATT(nid,id_SST ,"title","Temperature superficielle de la mer"),fnam)
252	CALL ncerr(NF90_PUT_ATT(nid,id_BILS,"title","Reference flux de chaleur au sol"),fnam)
253	CALL ncerr(NF90_PUT_ATT(nid,id_ALB, "title","Albedo a la surface"),fnam)
254	CALL ncerr(NF90_PUT_ATT(nid,id_RUG, "title","Rugosite"),fnam)
255
256	call nf95_put_att(nid, varid_longitude, "standard_name", "longitude")
257	call nf95_put_att(nid, varid_longitude, "units", "degrees_east")
258
259	call nf95_put_att(nid, varid_latitude, "standard_name", "latitude")
260	call nf95_put_att(nid, varid_latitude, "units", "degrees_north")
261
262	CALL ncerr(NF90_ENDDEF(nid),fnam)
263
264	!--- Variables saving
265	CALL ncerr(NF90_PUT_VAR(nid,id_tim,[(REAL(k),k=1,ndays)]),fnam)
266	CALL ncerr(NF90_PUT_VAR(nid,id_FOCE,pctsrf_t(:,is_oce,:),[1,1],[klon,ndays]),fnam)
267	CALL ncerr(NF90_PUT_VAR(nid,id_FSIC,pctsrf_t(:,is_sic,:),[1,1],[klon,ndays]),fnam)
268	CALL ncerr(NF90_PUT_VAR(nid,id_FTER,pctsrf_t(:,is_ter,:),[1,1],[klon,ndays]),fnam)
269	CALL ncerr(NF90_PUT_VAR(nid,id_FLIC,pctsrf_t(:,is_lic,:),[1,1],[klon,ndays]),fnam)
270	CALL ncerr(NF90_PUT_VAR(nid,id_SST ,phy_sst(:,:),[1,1],[klon,ndays]),fnam)
271	CALL ncerr(NF90_PUT_VAR(nid,id_BILS,phy_bil(:,:),[1,1],[klon,ndays]),fnam)
272	CALL ncerr(NF90_PUT_VAR(nid,id_ALB ,phy_alb(:,:),[1,1],[klon,ndays]),fnam)
273	CALL ncerr(NF90_PUT_VAR(nid,id_RUG ,phy_rug(:,:),[1,1],[klon,ndays]),fnam)
274	call nf95_put_var(nid, varid_longitude, longitude_deg)
275	call nf95_put_var(nid, varid_latitude, latitude_deg)
276
277	CALL ncerr(NF90_CLOSE(nid),fnam)
278
279	CALL msg(6,'Ecriture du fichier limit : fin')
280
281	DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug)
282
283
284	!===============================================================================
285	!
286	CONTAINS
287	!
288	!===============================================================================
289
290
291	!-------------------------------------------------------------------------------
292	!
293	SUBROUTINE get_2Dfield(fnam, varname, mode, ndays, champo, flag, mask)
294	!
295	!-----------------------------------------------------------------------------
296	! Comments:
297	! There are two assumptions concerning the NetCDF files, that are satisfied
298	! with files that are conforming NC convention:
299	! 1) The last dimension of the variables used is the time record.
300	! 2) Dimensional variables have the same names as corresponding dimensions.
301	!-----------------------------------------------------------------------------
302	USE netcdf, ONLY: NF90_OPEN, NF90_INQ_VARID, NF90_INQUIRE_VARIABLE, &
303	NF90_CLOSE, NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION, NF90_GET_VAR, &
304	NF90_GET_ATT
305	USE pchsp_95_m, only: pchsp_95
306	USE pchfe_95_m, only: pchfe_95
307	USE arth_m, only: arth
308	USE indice_sol_mod
309
310	IMPLICIT NONE
311	include "dimensions.h"
312	include "paramet.h"
313	include "comgeom2.h"
314	!-----------------------------------------------------------------------------
315	! Arguments:
316	CHARACTER(LEN=*), INTENT(IN) :: fnam ! NetCDF file name
317	CHARACTER(LEN=10), INTENT(IN) :: varname ! NetCDF variable name
318	CHARACTER(LEN=3), INTENT(IN) :: mode ! RUG, SIC, SST or ALB
319	INTEGER, INTENT(IN) :: ndays ! current year number of days
320	REAL, POINTER, DIMENSION(:, :) :: champo ! output field = f(t)
321	LOGICAL, OPTIONAL, INTENT(IN) :: flag ! extrapol. (SST) old ice (SIC)
322	REAL, OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask
323	!------------------------------------------------------------------------------
324	! Local variables:
325	!--- NetCDF
326	INTEGER :: ncid, varid ! NetCDF identifiers
327	CHARACTER(LEN=30) :: dnam ! dimension name
328	!--- dimensions
329	INTEGER :: dids(4) ! NetCDF dimensions identifiers
330	REAL, ALLOCATABLE :: dlon_ini(:) ! initial longitudes vector
331	REAL, ALLOCATABLE :: dlat_ini(:) ! initial latitudes vector
332	REAL, POINTER :: dlon(:), dlat(:) ! reordered lon/lat vectors
333	!--- fields
334	INTEGER :: imdep, jmdep, lmdep ! dimensions of 'champ'
335	REAL, ALLOCATABLE :: champ(:,:) ! wanted field on initial grid
336	REAL, ALLOCATABLE :: yder(:), timeyear(:)
337	REAL :: champint(iim,jjp1) ! interpolated field
338	REAL, ALLOCATABLE :: champtime(:,:,:)
339	REAL, ALLOCATABLE :: champan(:,:,:)
340	!--- input files
341	CHARACTER(LEN=20) :: fnam_m, fnam_p ! previous/next files names
342	CHARACTER(LEN=20) :: cal_in ! calendar
343	CHARACTER(LEN=20) :: unit_sic ! attribute unit in sea-ice file
344	INTEGER :: ndays_in ! number of days
345	!--- misc
346	INTEGER :: i, j, k, l, ll ! loop counters
347	REAL, ALLOCATABLE :: work(:,:) ! used for extrapolation
348	CHARACTER(LEN=128):: title, mess ! for messages
349	LOGICAL :: extrp ! flag for extrapolation
350	REAL :: chmin, chmax
351	INTEGER ierr, idx
352	integer n_extrap ! number of extrapolated points
353	logical skip
354
355	!------------------------------------------------------------------------------
356	!---Variables depending on keyword 'mode' -------------------------------------
357	NULLIFY(champo)
358
359	SELECT CASE(mode)
360	CASE('RUG'); title='Rugosite'
361	CASE('SIC'); title='Sea-ice'
362	CASE('SST'); title='SST'
363	CASE('ALB'); title='Albedo'
364	END SELECT
365	extrp=.FALSE.; IF(PRESENT(flag).AND.mode=='SST') extrp=flag
366	idx=INDEX(fnam,'.nc')-1
367
368	!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
369	CALL msg(5,' Now reading file : '//TRIM(fnam))
370	CALL ncerr(NF90_OPEN(fnam, NF90_NOWRITE, ncid),fnam)
371	CALL ncerr(NF90_INQ_VARID(ncid, trim(varname), varid),fnam)
372	CALL ncerr(NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids),fnam)
373
374	!--- Read unit for sea-ice variable only
375	IF (mode=='SIC') THEN
376	IF(NF90_GET_ATT(ncid, varid, 'units', unit_sic)/=NF90_NOERR) THEN
377	CALL msg(5,'No unit in sea-ice file. Take percentage as default value')
378	unit_sic='X'
379	ELSE
380	CALL msg(5,'Sea-ice cover has unit='//TRIM(unit_sic))
381	END IF
382	END IF
383
384	!--- Longitude
385	CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(1), name=dnam, len=imdep),fnam)
386	ALLOCATE(dlon_ini(imdep), dlon(imdep))
387	CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
388	CALL ncerr(NF90_GET_VAR(ncid, varid, dlon_ini), fnam)
389	CALL msg(5,'variable '//TRIM(dnam)//' dimension ', imdep)
390
391	!--- Latitude
392	CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(2), name=dnam, len=jmdep),fnam)
393	ALLOCATE(dlat_ini(jmdep), dlat(jmdep))
394	CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
395	CALL ncerr(NF90_GET_VAR(ncid, varid, dlat_ini), fnam)
396	CALL msg(5,'variable '//TRIM(dnam)//' dimension ', jmdep)
397
398	!--- Time (variable is not needed - it is rebuilt - but calendar is)
399	CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(3), name=dnam, len=lmdep), fnam)
400	ALLOCATE(timeyear(lmdep+2))
401	CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
402	cal_in=' '
403	IF(NF90_GET_ATT(ncid, varid, 'calendar', cal_in)/=NF90_NOERR) THEN
404	SELECT CASE(mode)
405	CASE('RUG', 'ALB'); cal_in='360d'
406	CASE('SIC', 'SST'); cal_in='gregorian'
407	END SELECT
408	CALL msg(5,'WARNING: missing "calendar" attribute for "time" in '&
409	&//TRIM(fnam)//'. Choosing default value.')
410	END IF
411	CALL msg(5,'var, calendar, dim: '//TRIM(dnam)//' '//TRIM(cal_in), lmdep)
412
413	!--- Determining input file number of days, depending on calendar
414	ndays_in=year_len(anneeref, cal_in)
415
416	!--- Rebuilding input time vector (field from input file might be unreliable)
417	IF(lmdep==12) THEN
418	timeyear=mid_month(anneeref, cal_in, ndays_in)
419	CALL msg(0,'Monthly input file(s) for '//TRIM(title)//'.')
420	ELSE IF(lmdep==ndays_in) THEN
421	timeyear=[(REAL(k)-0.5,k=0,ndays_in+1)]
422	CALL msg(0,'Daily input file (no time interpolation).')
423	ELSE
424	WRITE(mess,'(a,i3,a,i3,a)')'Mismatching input file: found',lmdep, &
425	' records, 12/',ndays_in,' (monthly/daily needed).'
426	CALL abort_physic('mid_months',TRIM(mess),1)
427	END IF
428
429	!--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
430	ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep+2))
431	IF(extrp) ALLOCATE(work(imdep, jmdep))
432	CALL msg(5,'')
433	CALL msg(5,'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, ' CHAMPS.')
434	CALL ncerr(NF90_INQ_VARID(ncid, varname, varid), fnam)
435	DO l=1, lmdep
436	CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,l],[imdep,jmdep,1]),fnam)
437	CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.)
438	IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
439	IF(l==1) THEN
440	CALL msg(5,"----------------------------------------------------------")
441	CALL msg(5,"$$$ Interpolation barycentrique pour "//TRIM(title)//" $$$")
442	CALL msg(5,"----------------------------------------------------------")
443	END IF
444	IF(mode=='RUG') champ=LOG(champ)
445	CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint)
446	IF(mode=='RUG') THEN
447	champint=EXP(champint)
448	WHERE(NINT(mask)/=1) champint=0.001
449	END IF
450	champtime(:, :, l+1)=champint
451	END DO
452	CALL ncerr(NF90_CLOSE(ncid), fnam)
453
454	!--- FIRST RECORD: LAST ONE OF PREVIOUS YEAR (CURRENT YEAR IF UNAVAILABLE)
455	fnam_m=fnam(1:idx)//'_m.nc'
456	IF(NF90_OPEN(fnam_m,NF90_NOWRITE,ncid)==NF90_NOERR) THEN
457	CALL msg(0,'Reading previous year file ("'//TRIM(fnam_m)//'") last record for '//TRIM(title))
458	CALL ncerr(NF90_INQ_VARID(ncid, varname, varid),fnam_m)
459	CALL ncerr(NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids),fnam_m)
460	CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(3), len=l), fnam_m)
461	CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,l],[imdep,jmdep,1]),fnam_m)
462	CALL ncerr(NF90_CLOSE(ncid), fnam_m)
463	CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.)
464	IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
465	IF(mode=='RUG') champ=LOG(champ)
466	CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint)
467	IF(mode=='RUG') THEN
468	champint=EXP(champint)
469	WHERE(NINT(mask)/=1) champint=0.001
470	END IF
471	champtime(:, :, 1)=champint
472	ELSE
473	CALL msg(0,'Using current year file ("'//TRIM(fnam)//'") last record for '//TRIM(title))
474	champtime(:, :, 1)=champtime(:, :, lmdep+1)
475	END IF
476
477	!--- LAST RECORD: FIRST ONE OF NEXT YEAR (CURRENT YEAR IF UNAVAILABLE)
478	fnam_p=fnam(1:idx)//'_p.nc'
479	IF(NF90_OPEN(fnam_p,NF90_NOWRITE,ncid)==NF90_NOERR) THEN
480	CALL msg(0,'Reading previous year file ("'//TRIM(fnam_p)//'") first record for '//TRIM(title))
481	CALL ncerr(NF90_INQ_VARID(ncid, varname, varid),fnam_p)
482	CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,1],[imdep,jmdep,1]),fnam_p)
483	CALL ncerr(NF90_CLOSE(ncid), fnam_p)
484	CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.)
485	IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
486	IF(mode=='RUG') champ=LOG(champ)
487	CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint)
488	IF(mode=='RUG') THEN
489	champint=EXP(champint)
490	WHERE(NINT(mask)/=1) champint=0.001
491	END IF
492	champtime(:, :, lmdep+2)=champint
493	ELSE
494	CALL msg(0,'Using current year file ("'//TRIM(fnam)//'") first record for '//TRIM(title))
495	champtime(:, :, lmdep+2)=champtime(:, :, 2)
496	END IF
497	DEALLOCATE(dlon_ini, dlat_ini, dlon, dlat, champ)
498	IF(extrp) DEALLOCATE(work)
499
500	!--- TIME INTERPOLATION ------------------------------------------------------
501	IF(prt_level>0) THEN
502	IF(ndays/=ndays_in) THEN
503	WRITE(lunout, *)
504	WRITE(lunout,*)'DIFFERENTES LONGEURS D ANNEES:'
505	WRITE(lunout,*)' Dans le fichier d entree: ',ndays_in
506	WRITE(lunout,*)' Dans le fichier de sortie: ',ndays
507	END IF
508	IF(lmdep==ndays_in) THEN
509	WRITE(lunout, *)
510	WRITE(lunout, *)'PAS D INTERPOLATION TEMPORELLE.'
511	WRITE(lunout, *)' Fichier journalier en entree.'
512	ELSE
513	WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
514	WRITE(lunout, *)' Vecteur temps en entree: ', timeyear
515	WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays-1
516	END IF
517	END IF
518	ALLOCATE(yder(lmdep+2), champan(iip1, jjp1, ndays))
519	IF(lmdep==ndays_in) THEN
520	champan(1:iim,:,:)=champtime
521	ELSE
522	skip = .false.
523	n_extrap = 0
524	DO j=1, jjp1
525	DO i=1, iim
526	yder = pchsp_95(timeyear, champtime(i, j, :), ibeg=2, iend=2, &
527	vc_beg=0., vc_end=0.)
528	CALL pchfe_95(timeyear, champtime(i, j, :), yder, skip, &
529	arth(0., real(ndays_in) / ndays, ndays), champan(i, j, :), ierr)
530	if (ierr < 0) stop 1
531	n_extrap = n_extrap + ierr
532	END DO
533	END DO
534	IF(n_extrap /= 0) WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
535	END IF
536	champan(iip1, :, :)=champan(1, :, :)
537	DEALLOCATE(yder, champtime, timeyear)
538
539	!--- Checking the result
540	DO j=1, jjp1
541	CALL minmax(iip1, champan(1, j, 10), chmin, chmax)
542	IF (prt_level>5) WRITE(lunout, *)' ',TRIM(title),' au temps 10 ', chmin, chmax, j
543	END DO
544
545	!--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
546	IF(mode=='SST') THEN
547	CALL msg(5,'Filtrage de la SST: SST >= 271.38')
548	WHERE(champan<271.38) champan=271.38
549	END IF
550
551	!--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
552	IF(mode=='SIC') THEN
553	CALL msg(5,'Filtrage de la SIC: 0.0 < Sea-ice < 1.0')
554
555	IF (unit_sic=='1') THEN
556	! Nothing to be done for sea-ice field is already in fraction of 1
557	! This is the case for sea-ice in file cpl_atm_sic.nc
558	CALL msg(5,'Sea-ice field already in fraction of 1')
559	ELSE
560	! Convert sea ice from percentage to fraction of 1
561	CALL msg(5,'Transformt sea-ice field from percentage to fraction of 1.')
562	champan(:, :, :)=champan(:, :, :)/100.
563	END IF
564
565	champan(iip1, :, :)=champan(1, :, :)
566	WHERE(champan>1.0) champan=1.0
567	WHERE(champan<0.0) champan=0.0
568	END IF
569
570	!--- DYNAMICAL TO PHYSICAL GRID ----------------------------------------------
571	ALLOCATE(champo(klon, ndays))
572	DO k=1, ndays
573	CALL gr_dyn_fi(1, iip1, jjp1, klon, champan(1, 1, k), champo(1, k))
574	END DO
575	DEALLOCATE(champan)
576
577	END SUBROUTINE get_2Dfield
578	!
579	!-------------------------------------------------------------------------------
580
581
582	!-------------------------------------------------------------------------------
583	!
584	SUBROUTINE start_init_orog0(lon_in,lat_in,phis,masque)
585	!
586	!-------------------------------------------------------------------------------
587	USE grid_noro_m, ONLY: grid_noro0
588	IMPLICIT NONE
589	!===============================================================================
590	! Purpose: Compute "phis" just like it would be in start_init_orog.
591	!===============================================================================
592	! Arguments:
593	REAL, INTENT(IN) :: lon_in(:), lat_in(:) ! dim (iml) (jml)
594	REAL, INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml)
595	!-------------------------------------------------------------------------------
596	! Local variables:
597	CHARACTER(LEN=256) :: modname="start_init_orog0"
598	INTEGER :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1)
599	REAL :: lev(1), date, dt, deg2rad
600	REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:)
601	REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:), lat_rel(:,:)
602	!-------------------------------------------------------------------------------
603	iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml")
604	jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml")
605	IF(iml/=iip1) CALL abort_gcm(TRIM(modname),'iml/=iip1',1)
606	IF(jml/=jjp1) CALL abort_gcm(TRIM(modname),'jml/=jjp1',1)
607	pi=2.0*ASIN(1.0); deg2rad=pi/180.0
608	IF(ANY(phis/=-99999.)) RETURN !--- phis ALREADY KNOWN
609
610	!--- HIGH RESOLUTION OROGRAPHY
611	CALL flininfo(frelf, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
612
613	ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
614	CALL flinopen(frelf, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel, &
615	lev, ttm_tmp, itau, date, dt, fid)
616	ALLOCATE(relief_hi(iml_rel,jml_rel))
617	CALL flinget(fid, vrel, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
618	CALL flinclo(fid)
619
620	!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
621	ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
622	lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
623	lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
624
625	!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
626	ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
627	CALL conf_dat2d(vrel, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.)
628	DEALLOCATE(lon_ini,lat_ini)
629
630	!--- COMPUTING SURFACE GEOPOTENTIAL USING ROUTINE grid_noro0
631	WRITE(lunout,*)
632	WRITE(lunout,)' Compute surface geopotential *'
633
634	!--- CALL OROGRAPHY MODULE (REDUCED VERSION) TO COMPUTE FIELDS
635	CALL grid_noro0(lon_rad, lat_rad, relief_hi, lon_in, lat_in, phis, masque)
636	phis = phis * 9.81
637	phis(iml,:) = phis(1,:)
638	DEALLOCATE(relief_hi,lon_rad,lat_rad)
639
640	END SUBROUTINE start_init_orog0
641	!
642	!-------------------------------------------------------------------------------
643
644
645	!-------------------------------------------------------------------------------
646	!
647	FUNCTION year_len(y,cal_in)
648	!
649	!-------------------------------------------------------------------------------
650	IMPLICIT NONE
651	!-------------------------------------------------------------------------------
652	! Arguments:
653	INTEGER :: year_len
654	INTEGER, INTENT(IN) :: y
655	CHARACTER(LEN=*), INTENT(IN) :: cal_in
656	!-------------------------------------------------------------------------------
657	! Local variables:
658	CHARACTER(LEN=20) :: cal_out ! calendar (for outputs)
659	!-------------------------------------------------------------------------------
660	!--- Getting the input calendar to reset at the end of the function
661	CALL ioget_calendar(cal_out)
662
663	!--- Unlocking calendar and setting it to wanted one
664	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
665
666	!--- Getting the number of days in this year
667	year_len=ioget_year_len(y)
668
669	!--- Back to original calendar
670	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
671
672	END FUNCTION year_len
673	!
674	!-------------------------------------------------------------------------------
675
676
677	!-------------------------------------------------------------------------------
678	!
679	FUNCTION mid_month(y,cal_in,ndays_out)
680	!
681	!-------------------------------------------------------------------------------
682	IMPLICIT NONE
683	!-------------------------------------------------------------------------------
684	! Arguments:
685	INTEGER, INTENT(IN) :: y ! year
686	CHARACTER(LEN=*), INTENT(IN) :: cal_in ! calendar
687	INTEGER, INTENT(IN) :: ndays_out ! days number
688	REAL, DIMENSION(14) :: mid_month ! mid-bins times
689	!-------------------------------------------------------------------------------
690	! Local variables:
691	CHARACTER(LEN=99) :: mess ! error message
692	CHARACTER(LEN=20) :: cal_out ! output calendar
693	INTEGER, DIMENSION(14) :: tlen ! months lengths (days)
694	INTEGER :: m ! months counter
695	INTEGER :: nd ! number of days
696	!-------------------------------------------------------------------------------
697	!--- Save the input calendar to reset it at the end of the function
698	CALL ioget_calendar(cal_out)
699
700	!--- Unlock calendar and set it to "cal_in"
701	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
702
703	!--- Get the length of each month
704	tlen(1 )=ioget_mon_len(y-1,12)
705	DO m=1,12; tlen(m+1)=ioget_mon_len(y,m); END DO
706	tlen(14)=ioget_mon_len(y+1, 1)
707
708	!--- Mid-bins times
709	mid_month(1)=-0.5*REAL(tlen(1))
710	DO m=2,14; mid_month(m)=mid_month(m-1)+0.5*REAL(tlen(m-1)+tlen(m)); END DO
711
712	!--- Back to original calendar
713	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
714
715	END FUNCTION mid_month
716	!
717	!-------------------------------------------------------------------------------
718
719
720
721	!-------------------------------------------------------------------------------
722	!
723	SUBROUTINE msg(lev,str1,i,str2)
724	!
725	!-------------------------------------------------------------------------------
726	! Arguments:
727	INTEGER, INTENT(IN) :: lev
728	CHARACTER(LEN=*), INTENT(IN) :: str1
729	INTEGER, OPTIONAL, INTENT(IN) :: i
730	CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: str2
731	!-------------------------------------------------------------------------------
732	IF(prt_level>lev) THEN
733	IF(PRESENT(str2)) THEN
734	WRITE(lunout,*) TRIM(str1), i, TRIM(str2)
735	ELSE IF(PRESENT(i)) THEN
736	WRITE(lunout,*) TRIM(str1), i
737	ELSE
738	WRITE(lunout,*) TRIM(str1)
739	END IF
740	END IF
741
742	END SUBROUTINE msg
743	!
744	!-------------------------------------------------------------------------------
745
746
747	!-------------------------------------------------------------------------------
748	!
749	SUBROUTINE ncerr(ncres,fnam)
750	!
751	!-------------------------------------------------------------------------------
752	! Purpose: NetCDF errors handling.
753	!-------------------------------------------------------------------------------
754	USE netcdf, ONLY : NF90_NOERR, NF90_STRERROR
755	IMPLICIT NONE
756	!-------------------------------------------------------------------------------
757	! Arguments:
758	INTEGER, INTENT(IN) :: ncres
759	CHARACTER(LEN=*), INTENT(IN) :: fnam
760	!-------------------------------------------------------------------------------
761	IF(ncres/=NF90_NOERR) THEN
762	WRITE(lunout,*)'Problem with file '//TRIM(fnam)//' in routine limit_netcdf.'
763	CALL abort_physic('limit_netcdf',NF90_STRERROR(ncres),1)
764	END IF
765
766	END SUBROUTINE ncerr
767	!
768	!-------------------------------------------------------------------------------
769
770	#endif
771	! of #ifndef CPP_1D
772	END SUBROUTINE limit_netcdf
773
774	END MODULE limit
775	!
776	!*******************************************************************************
777

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

Download in other formats:

Original Format