Context Navigation

limit_netcdf.F90 @ 2665

Last change on this file since 2665 was 2665, checked in by dcugnet, 8 years ago
A (re)startphy.nc file (standard name: "startphy0.nc") can be read by ce0l to get land mask, so mask can be defined (in decreasing priority order) from: 1) "o2a.nc file" if this file is found 2) "startphy0.nc" if this file is found 3) "Relief.nc" otherwise Sub-cell scales parameters for orographic gravity waves can be read from file "oro_params.nc" if the configuration key "read_orop" is TRUE. The effect is to bypass the "grid_noro" routine in ce0l, so that any pre-defined mask (from o2a.nc or startphy0.nc) is then overwritten. The gcm stops if the "limit.nc" records number differs from the current year number of days. A warning is issued in case the gcm calendar does not match the time axis attribute "calendar" (if available) from the "limit.nc" file. This attribute is now added to the "limit.nc" time axis. Few simplifications in grid_noro Few parameters changes in acama_gwd and flott_gwd. Variable d_u can be saved in the outputs.
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: 33.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	USE comconst_mod, ONLY: pi
73	USE phys_cal_mod, ONLY: calend
74	IMPLICIT NONE
75	!-------------------------------------------------------------------------------
76	! Arguments:
77	include "iniprint.h"
78	include "dimensions.h"
79	include "paramet.h"
80	REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: masque ! land mask
81	REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: phis ! ground geopotential
82	LOGICAL, INTENT(IN) :: extrap ! SST extrapolation flag
83	!-------------------------------------------------------------------------------
84	! Local variables:
85	include "comgeom2.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_tim, "calendar",calend),fnam)
248	CALL ncerr(NF90_PUT_ATT(nid,id_FOCE,"title","Fraction ocean"),fnam)
249	CALL ncerr(NF90_PUT_ATT(nid,id_FSIC,"title","Fraction glace de mer"),fnam)
250	CALL ncerr(NF90_PUT_ATT(nid,id_FTER,"title","Fraction terre"),fnam)
251	CALL ncerr(NF90_PUT_ATT(nid,id_FLIC,"title","Fraction land ice"),fnam)
252	CALL ncerr(NF90_PUT_ATT(nid,id_SST ,"title","Temperature superficielle de la mer"),fnam)
253	CALL ncerr(NF90_PUT_ATT(nid,id_BILS,"title","Reference flux de chaleur au sol"),fnam)
254	CALL ncerr(NF90_PUT_ATT(nid,id_ALB, "title","Albedo a la surface"),fnam)
255	CALL ncerr(NF90_PUT_ATT(nid,id_RUG, "title","Rugosite"),fnam)
256
257	call nf95_put_att(nid, varid_longitude, "standard_name", "longitude")
258	call nf95_put_att(nid, varid_longitude, "units", "degrees_east")
259
260	call nf95_put_att(nid, varid_latitude, "standard_name", "latitude")
261	call nf95_put_att(nid, varid_latitude, "units", "degrees_north")
262
263	CALL ncerr(NF90_ENDDEF(nid),fnam)
264
265	!--- Variables saving
266	CALL ncerr(NF90_PUT_VAR(nid,id_tim,[(REAL(k),k=1,ndays)]),fnam)
267	CALL ncerr(NF90_PUT_VAR(nid,id_FOCE,pctsrf_t(:,is_oce,:),[1,1],[klon,ndays]),fnam)
268	CALL ncerr(NF90_PUT_VAR(nid,id_FSIC,pctsrf_t(:,is_sic,:),[1,1],[klon,ndays]),fnam)
269	CALL ncerr(NF90_PUT_VAR(nid,id_FTER,pctsrf_t(:,is_ter,:),[1,1],[klon,ndays]),fnam)
270	CALL ncerr(NF90_PUT_VAR(nid,id_FLIC,pctsrf_t(:,is_lic,:),[1,1],[klon,ndays]),fnam)
271	CALL ncerr(NF90_PUT_VAR(nid,id_SST ,phy_sst(:,:),[1,1],[klon,ndays]),fnam)
272	CALL ncerr(NF90_PUT_VAR(nid,id_BILS,phy_bil(:,:),[1,1],[klon,ndays]),fnam)
273	CALL ncerr(NF90_PUT_VAR(nid,id_ALB ,phy_alb(:,:),[1,1],[klon,ndays]),fnam)
274	CALL ncerr(NF90_PUT_VAR(nid,id_RUG ,phy_rug(:,:),[1,1],[klon,ndays]),fnam)
275	call nf95_put_var(nid, varid_longitude, longitude_deg)
276	call nf95_put_var(nid, varid_latitude, latitude_deg)
277
278	CALL ncerr(NF90_CLOSE(nid),fnam)
279
280	CALL msg(6,'Ecriture du fichier limit : fin')
281
282	DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug)
283
284
285	!===============================================================================
286	!
287	CONTAINS
288	!
289	!===============================================================================
290
291
292	!-------------------------------------------------------------------------------
293	!
294	SUBROUTINE get_2Dfield(fnam, varname, mode, ndays, champo, flag, mask)
295	!
296	!-----------------------------------------------------------------------------
297	! Comments:
298	! There are two assumptions concerning the NetCDF files, that are satisfied
299	! with files that are conforming NC convention:
300	! 1) The last dimension of the variables used is the time record.
301	! 2) Dimensional variables have the same names as corresponding dimensions.
302	!-----------------------------------------------------------------------------
303	USE netcdf, ONLY: NF90_OPEN, NF90_INQ_VARID, NF90_INQUIRE_VARIABLE, &
304	NF90_CLOSE, NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION, NF90_GET_VAR, &
305	NF90_GET_ATT
306	USE pchsp_95_m, only: pchsp_95
307	USE pchfe_95_m, only: pchfe_95
308	USE arth_m, only: arth
309	USE indice_sol_mod
310
311	IMPLICIT NONE
312	include "dimensions.h"
313	include "paramet.h"
314	include "comgeom2.h"
315	!-----------------------------------------------------------------------------
316	! Arguments:
317	CHARACTER(LEN=*), INTENT(IN) :: fnam ! NetCDF file name
318	CHARACTER(LEN=10), INTENT(IN) :: varname ! NetCDF variable name
319	CHARACTER(LEN=3), INTENT(IN) :: mode ! RUG, SIC, SST or ALB
320	INTEGER, INTENT(IN) :: ndays ! current year number of days
321	REAL, POINTER, DIMENSION(:, :) :: champo ! output field = f(t)
322	LOGICAL, OPTIONAL, INTENT(IN) :: flag ! extrapol. (SST) old ice (SIC)
323	REAL, OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask
324	!------------------------------------------------------------------------------
325	! Local variables:
326	!--- NetCDF
327	INTEGER :: ncid, varid ! NetCDF identifiers
328	CHARACTER(LEN=30) :: dnam ! dimension name
329	!--- dimensions
330	INTEGER :: dids(4) ! NetCDF dimensions identifiers
331	REAL, ALLOCATABLE :: dlon_ini(:) ! initial longitudes vector
332	REAL, ALLOCATABLE :: dlat_ini(:) ! initial latitudes vector
333	REAL, POINTER :: dlon(:), dlat(:) ! reordered lon/lat vectors
334	!--- fields
335	INTEGER :: imdep, jmdep, lmdep ! dimensions of 'champ'
336	REAL, ALLOCATABLE :: champ(:,:) ! wanted field on initial grid
337	REAL, ALLOCATABLE :: yder(:), timeyear(:)
338	REAL :: champint(iim,jjp1) ! interpolated field
339	REAL, ALLOCATABLE :: champtime(:,:,:)
340	REAL, ALLOCATABLE :: champan(:,:,:)
341	!--- input files
342	CHARACTER(LEN=20) :: fnam_m, fnam_p ! previous/next files names
343	CHARACTER(LEN=20) :: cal_in ! calendar
344	CHARACTER(LEN=20) :: unit_sic ! attribute unit in sea-ice file
345	INTEGER :: ndays_in ! number of days
346	!--- misc
347	INTEGER :: i, j, k, l, ll ! loop counters
348	REAL, ALLOCATABLE :: work(:,:) ! used for extrapolation
349	CHARACTER(LEN=128):: title, mess ! for messages
350	LOGICAL :: extrp ! flag for extrapolation
351	REAL :: chmin, chmax
352	INTEGER ierr, idx
353	integer n_extrap ! number of extrapolated points
354	logical skip
355
356	!------------------------------------------------------------------------------
357	!---Variables depending on keyword 'mode' -------------------------------------
358	NULLIFY(champo)
359
360	SELECT CASE(mode)
361	CASE('RUG'); title='Rugosite'
362	CASE('SIC'); title='Sea-ice'
363	CASE('SST'); title='SST'
364	CASE('ALB'); title='Albedo'
365	END SELECT
366	extrp=.FALSE.; IF(PRESENT(flag).AND.mode=='SST') extrp=flag
367	idx=INDEX(fnam,'.nc')-1
368
369	!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
370	CALL msg(5,' Now reading file : '//TRIM(fnam))
371	CALL ncerr(NF90_OPEN(fnam, NF90_NOWRITE, ncid),fnam)
372	CALL ncerr(NF90_INQ_VARID(ncid, trim(varname), varid),fnam)
373	CALL ncerr(NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids),fnam)
374
375	!--- Read unit for sea-ice variable only
376	IF (mode=='SIC') THEN
377	IF(NF90_GET_ATT(ncid, varid, 'units', unit_sic)/=NF90_NOERR) THEN
378	CALL msg(5,'No unit in sea-ice file. Take percentage as default value')
379	unit_sic='X'
380	ELSE
381	CALL msg(5,'Sea-ice cover has unit='//TRIM(unit_sic))
382	END IF
383	END IF
384
385	!--- Longitude
386	CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(1), name=dnam, len=imdep),fnam)
387	ALLOCATE(dlon_ini(imdep), dlon(imdep))
388	CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
389	CALL ncerr(NF90_GET_VAR(ncid, varid, dlon_ini), fnam)
390	CALL msg(5,'variable '//TRIM(dnam)//' dimension ', imdep)
391
392	!--- Latitude
393	CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(2), name=dnam, len=jmdep),fnam)
394	ALLOCATE(dlat_ini(jmdep), dlat(jmdep))
395	CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
396	CALL ncerr(NF90_GET_VAR(ncid, varid, dlat_ini), fnam)
397	CALL msg(5,'variable '//TRIM(dnam)//' dimension ', jmdep)
398
399	!--- Time (variable is not needed - it is rebuilt - but calendar is)
400	CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(3), name=dnam, len=lmdep), fnam)
401	ALLOCATE(timeyear(lmdep+2))
402	CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
403	cal_in=' '
404	IF(NF90_GET_ATT(ncid, varid, 'calendar', cal_in)/=NF90_NOERR) THEN
405	SELECT CASE(mode)
406	CASE('RUG', 'ALB'); cal_in='360d'
407	CASE('SIC', 'SST'); cal_in='gregorian'
408	END SELECT
409	CALL msg(5,'WARNING: missing "calendar" attribute for "time" in '&
410	&//TRIM(fnam)//'. Choosing default value.')
411	END IF
412	CALL strclean(cal_in) !--- REMOVE (WEIRD) NULL CHARACTERS
413	CALL msg(5,'var, calendar, dim: '//TRIM(dnam)//' '//TRIM(cal_in), lmdep)
414
415	!--- Determining input file number of days, depending on calendar
416	ndays_in=year_len(anneeref, cal_in)
417
418	!--- Rebuilding input time vector (field from input file might be unreliable)
419	IF(lmdep==12) THEN
420	timeyear=mid_month(anneeref, cal_in, ndays_in)
421	CALL msg(0,'Monthly input file(s) for '//TRIM(title)//'.')
422	ELSE IF(lmdep==ndays_in) THEN
423	timeyear=[(REAL(k)-0.5,k=0,ndays_in+1)]
424	CALL msg(0,'Daily input file (no time interpolation).')
425	ELSE
426	WRITE(mess,'(a,i3,a,i3,a)')'Mismatching input file: found',lmdep, &
427	' records, 12/',ndays_in,' (monthly/daily needed).'
428	CALL abort_physic('mid_months',TRIM(mess),1)
429	END IF
430
431	!--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
432	ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep+2))
433	IF(extrp) ALLOCATE(work(imdep, jmdep))
434	CALL msg(5,'')
435	CALL msg(5,'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, ' CHAMPS.')
436	CALL ncerr(NF90_INQ_VARID(ncid, varname, varid), fnam)
437	DO l=1, lmdep
438	CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,l],[imdep,jmdep,1]),fnam)
439	CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.)
440	IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
441	IF(l==1) THEN
442	CALL msg(5,"----------------------------------------------------------")
443	CALL msg(5,"$$$ Interpolation barycentrique pour "//TRIM(title)//" $$$")
444	CALL msg(5,"----------------------------------------------------------")
445	END IF
446	IF(mode=='RUG') champ=LOG(champ)
447	CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint)
448	IF(mode=='RUG') THEN
449	champint=EXP(champint)
450	WHERE(NINT(mask)/=1) champint=0.001
451	END IF
452	champtime(:, :, l+1)=champint
453	END DO
454	CALL ncerr(NF90_CLOSE(ncid), fnam)
455
456	!--- FIRST RECORD: LAST ONE OF PREVIOUS YEAR (CURRENT YEAR IF UNAVAILABLE)
457	fnam_m=fnam(1:idx)//'_m.nc'
458	IF(NF90_OPEN(fnam_m,NF90_NOWRITE,ncid)==NF90_NOERR) THEN
459	CALL msg(0,'Reading previous year file ("'//TRIM(fnam_m)//'") last record for '//TRIM(title))
460	CALL ncerr(NF90_INQ_VARID(ncid, varname, varid),fnam_m)
461	CALL ncerr(NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids),fnam_m)
462	CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(3), len=l), fnam_m)
463	CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,l],[imdep,jmdep,1]),fnam_m)
464	CALL ncerr(NF90_CLOSE(ncid), fnam_m)
465	CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.)
466	IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
467	IF(mode=='RUG') champ=LOG(champ)
468	CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint)
469	IF(mode=='RUG') THEN
470	champint=EXP(champint)
471	WHERE(NINT(mask)/=1) champint=0.001
472	END IF
473	champtime(:, :, 1)=champint
474	ELSE
475	CALL msg(0,'Using current year file ("'//TRIM(fnam)//'") last record for '//TRIM(title))
476	champtime(:, :, 1)=champtime(:, :, lmdep+1)
477	END IF
478
479	!--- LAST RECORD: FIRST ONE OF NEXT YEAR (CURRENT YEAR IF UNAVAILABLE)
480	fnam_p=fnam(1:idx)//'_p.nc'
481	IF(NF90_OPEN(fnam_p,NF90_NOWRITE,ncid)==NF90_NOERR) THEN
482	CALL msg(0,'Reading next year file ("'//TRIM(fnam_p)//'") first record for '//TRIM(title))
483	CALL ncerr(NF90_INQ_VARID(ncid, varname, varid),fnam_p)
484	CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,1],[imdep,jmdep,1]),fnam_p)
485	CALL ncerr(NF90_CLOSE(ncid), fnam_p)
486	CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.)
487	IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
488	IF(mode=='RUG') champ=LOG(champ)
489	CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint)
490	IF(mode=='RUG') THEN
491	champint=EXP(champint)
492	WHERE(NINT(mask)/=1) champint=0.001
493	END IF
494	champtime(:, :, lmdep+2)=champint
495	ELSE
496	CALL msg(0,'Using current year file ("'//TRIM(fnam)//'") first record for '//TRIM(title))
497	champtime(:, :, lmdep+2)=champtime(:, :, 2)
498	END IF
499	DEALLOCATE(dlon_ini, dlat_ini, dlon, dlat, champ)
500	IF(extrp) DEALLOCATE(work)
501
502	!--- TIME INTERPOLATION ------------------------------------------------------
503	IF(prt_level>0) THEN
504	IF(ndays/=ndays_in) THEN
505	WRITE(lunout, *)
506	WRITE(lunout,*)'DIFFERENTES LONGEURS D ANNEES:'
507	WRITE(lunout,*)' Dans le fichier d entree: ',ndays_in
508	WRITE(lunout,*)' Dans le fichier de sortie: ',ndays
509	END IF
510	IF(lmdep==ndays_in) THEN
511	WRITE(lunout, *)
512	WRITE(lunout, *)'PAS D INTERPOLATION TEMPORELLE.'
513	WRITE(lunout, *)' Fichier journalier en entree.'
514	ELSE
515	WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
516	WRITE(lunout, *)' Vecteur temps en entree: ', timeyear
517	WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays-1
518	END IF
519	END IF
520	ALLOCATE(yder(lmdep+2), champan(iip1, jjp1, ndays))
521	IF(lmdep==ndays_in) THEN
522	champan(1:iim,:,:)=champtime
523	ELSE
524	skip = .false.
525	n_extrap = 0
526	DO j=1, jjp1
527	DO i=1, iim
528	yder = pchsp_95(timeyear, champtime(i, j, :), ibeg=2, iend=2, &
529	vc_beg=0., vc_end=0.)
530	CALL pchfe_95(timeyear, champtime(i, j, :), yder, skip, &
531	arth(0., real(ndays_in) / ndays, ndays), champan(i, j, :), ierr)
532	if (ierr < 0) stop 1
533	n_extrap = n_extrap + ierr
534	END DO
535	END DO
536	IF(n_extrap /= 0) WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
537	END IF
538	champan(iip1, :, :)=champan(1, :, :)
539	DEALLOCATE(yder, champtime, timeyear)
540
541	!--- Checking the result
542	DO j=1, jjp1
543	CALL minmax(iip1, champan(1, j, 10), chmin, chmax)
544	IF (prt_level>5) WRITE(lunout, *)' ',TRIM(title),' au temps 10 ', chmin, chmax, j
545	END DO
546
547	!--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
548	IF(mode=='SST') THEN
549	CALL msg(5,'Filtrage de la SST: SST >= 271.38')
550	WHERE(champan<271.38) champan=271.38
551	END IF
552
553	!--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
554	IF(mode=='SIC') THEN
555	CALL msg(5,'Filtrage de la SIC: 0.0 < Sea-ice < 1.0')
556
557	IF (unit_sic=='1') THEN
558	! Nothing to be done for sea-ice field is already in fraction of 1
559	! This is the case for sea-ice in file cpl_atm_sic.nc
560	CALL msg(5,'Sea-ice field already in fraction of 1')
561	ELSE
562	! Convert sea ice from percentage to fraction of 1
563	CALL msg(5,'Transformt sea-ice field from percentage to fraction of 1.')
564	champan(:, :, :)=champan(:, :, :)/100.
565	END IF
566
567	champan(iip1, :, :)=champan(1, :, :)
568	WHERE(champan>1.0) champan=1.0
569	WHERE(champan<0.0) champan=0.0
570	END IF
571
572	!--- DYNAMICAL TO PHYSICAL GRID ----------------------------------------------
573	ALLOCATE(champo(klon, ndays))
574	DO k=1, ndays
575	CALL gr_dyn_fi(1, iip1, jjp1, klon, champan(1, 1, k), champo(1, k))
576	END DO
577	DEALLOCATE(champan)
578
579	END SUBROUTINE get_2Dfield
580	!
581	!-------------------------------------------------------------------------------
582
583
584	!-------------------------------------------------------------------------------
585	!
586	SUBROUTINE start_init_orog0(lon_in,lat_in,phis,masque)
587	!
588	!-------------------------------------------------------------------------------
589	USE grid_noro_m, ONLY: grid_noro0
590	IMPLICIT NONE
591	!===============================================================================
592	! Purpose: Compute "phis" just like it would be in start_init_orog.
593	!===============================================================================
594	! Arguments:
595	REAL, INTENT(IN) :: lon_in(:), lat_in(:) ! dim (iml) (jml)
596	REAL, INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml)
597	!-------------------------------------------------------------------------------
598	! Local variables:
599	CHARACTER(LEN=256) :: modname="start_init_orog0"
600	INTEGER :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1)
601	REAL :: lev(1), date, dt, deg2rad
602	REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:)
603	REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:), lat_rel(:,:)
604	!-------------------------------------------------------------------------------
605	iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml")
606	jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml")
607	IF(iml/=iip1) CALL abort_gcm(TRIM(modname),'iml/=iip1',1)
608	IF(jml/=jjp1) CALL abort_gcm(TRIM(modname),'jml/=jjp1',1)
609	pi=2.0*ASIN(1.0); deg2rad=pi/180.0
610	IF(ANY(phis/=-99999.)) RETURN !--- phis ALREADY KNOWN
611
612	!--- HIGH RESOLUTION OROGRAPHY
613	CALL flininfo(frelf, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
614
615	ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
616	CALL flinopen(frelf, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel, &
617	lev, ttm_tmp, itau, date, dt, fid)
618	ALLOCATE(relief_hi(iml_rel,jml_rel))
619	CALL flinget(fid, vrel, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
620	CALL flinclo(fid)
621
622	!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
623	ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
624	lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
625	lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
626
627	!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
628	ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
629	CALL conf_dat2d(vrel, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.)
630	DEALLOCATE(lon_ini,lat_ini)
631
632	!--- COMPUTING SURFACE GEOPOTENTIAL USING ROUTINE grid_noro0
633	WRITE(lunout,*)
634	WRITE(lunout,)' Compute surface geopotential *'
635
636	!--- CALL OROGRAPHY MODULE (REDUCED VERSION) TO COMPUTE FIELDS
637	CALL grid_noro0(lon_rad, lat_rad, relief_hi, lon_in, lat_in, phis, masque)
638	phis = phis * 9.81
639	phis(iml,:) = phis(1,:)
640	DEALLOCATE(relief_hi,lon_rad,lat_rad)
641
642	END SUBROUTINE start_init_orog0
643	!
644	!-------------------------------------------------------------------------------
645
646
647	!-------------------------------------------------------------------------------
648	!
649	FUNCTION year_len(y,cal_in)
650	!
651	!-------------------------------------------------------------------------------
652	IMPLICIT NONE
653	!-------------------------------------------------------------------------------
654	! Arguments:
655	INTEGER :: year_len
656	INTEGER, INTENT(IN) :: y
657	CHARACTER(LEN=*), INTENT(IN) :: cal_in
658	!-------------------------------------------------------------------------------
659	! Local variables:
660	CHARACTER(LEN=20) :: cal_out ! calendar (for outputs)
661	!-------------------------------------------------------------------------------
662	!--- Getting the input calendar to reset at the end of the function
663	CALL ioget_calendar(cal_out)
664
665	!--- Unlocking calendar and setting it to wanted one
666	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
667
668	!--- Getting the number of days in this year
669	year_len=ioget_year_len(y)
670
671	!--- Back to original calendar
672	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
673
674	END FUNCTION year_len
675	!
676	!-------------------------------------------------------------------------------
677
678
679	!-------------------------------------------------------------------------------
680	!
681	FUNCTION mid_month(y,cal_in,ndays_out)
682	!
683	!-------------------------------------------------------------------------------
684	IMPLICIT NONE
685	!-------------------------------------------------------------------------------
686	! Arguments:
687	INTEGER, INTENT(IN) :: y ! year
688	CHARACTER(LEN=*), INTENT(IN) :: cal_in ! calendar
689	INTEGER, INTENT(IN) :: ndays_out ! days number
690	REAL, DIMENSION(14) :: mid_month ! mid-bins times
691	!-------------------------------------------------------------------------------
692	! Local variables:
693	CHARACTER(LEN=99) :: mess ! error message
694	CHARACTER(LEN=20) :: cal_out ! output calendar
695	INTEGER, DIMENSION(14) :: tlen ! months lengths (days)
696	INTEGER :: m ! months counter
697	INTEGER :: nd ! number of days
698	!-------------------------------------------------------------------------------
699	!--- Save the input calendar to reset it at the end of the function
700	CALL ioget_calendar(cal_out)
701
702	!--- Unlock calendar and set it to "cal_in"
703	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
704
705	!--- Get the length of each month
706	tlen(1 )=ioget_mon_len(y-1,12)
707	DO m=1,12; tlen(m+1)=ioget_mon_len(y,m); END DO
708	tlen(14)=ioget_mon_len(y+1, 1)
709
710	!--- Mid-bins times
711	mid_month(1)=-0.5*REAL(tlen(1))
712	DO m=2,14; mid_month(m)=mid_month(m-1)+0.5*REAL(tlen(m-1)+tlen(m)); END DO
713
714	!--- Back to original calendar
715	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
716
717	END FUNCTION mid_month
718	!
719	!-------------------------------------------------------------------------------
720
721
722
723	!-------------------------------------------------------------------------------
724	!
725	SUBROUTINE msg(lev,str1,i,str2)
726	!
727	!-------------------------------------------------------------------------------
728	! Arguments:
729	INTEGER, INTENT(IN) :: lev
730	CHARACTER(LEN=*), INTENT(IN) :: str1
731	INTEGER, OPTIONAL, INTENT(IN) :: i
732	CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: str2
733	!-------------------------------------------------------------------------------
734	IF(prt_level>lev) THEN
735	IF(PRESENT(str2)) THEN
736	WRITE(lunout,*) TRIM(str1), i, TRIM(str2)
737	ELSE IF(PRESENT(i)) THEN
738	WRITE(lunout,*) TRIM(str1), i
739	ELSE
740	WRITE(lunout,*) TRIM(str1)
741	END IF
742	END IF
743
744	END SUBROUTINE msg
745	!
746	!-------------------------------------------------------------------------------
747
748
749	!-------------------------------------------------------------------------------
750	!
751	SUBROUTINE ncerr(ncres,fnam)
752	!
753	!-------------------------------------------------------------------------------
754	! Purpose: NetCDF errors handling.
755	!-------------------------------------------------------------------------------
756	USE netcdf, ONLY : NF90_NOERR, NF90_STRERROR
757	IMPLICIT NONE
758	!-------------------------------------------------------------------------------
759	! Arguments:
760	INTEGER, INTENT(IN) :: ncres
761	CHARACTER(LEN=*), INTENT(IN) :: fnam
762	!-------------------------------------------------------------------------------
763	IF(ncres/=NF90_NOERR) THEN
764	WRITE(lunout,*)'Problem with file '//TRIM(fnam)//' in routine limit_netcdf.'
765	CALL abort_physic('limit_netcdf',NF90_STRERROR(ncres),1)
766	END IF
767
768	END SUBROUTINE ncerr
769	!
770	!-------------------------------------------------------------------------------
771
772
773	!-------------------------------------------------------------------------------
774	!
775	SUBROUTINE strclean(s)
776	!
777	!-------------------------------------------------------------------------------
778	IMPLICIT NONE
779	!-------------------------------------------------------------------------------
780	! Purpose: Remove tail null characters from the input string.
781	!-------------------------------------------------------------------------------
782	! Parameters:
783	CHARACTER(LEN=*), INTENT(INOUT) :: s
784	!-------------------------------------------------------------------------------
785	! Local variable:
786	INTEGER :: k
787	!-------------------------------------------------------------------------------
788	k=LEN_TRIM(s); DO WHILE(ICHAR(s(k:k))==0); s(k:k)=' '; k=LEN_TRIM(s); END DO
789
790	END SUBROUTINE strclean
791	!
792	!-------------------------------------------------------------------------------
793
794	#endif
795	! of #ifndef CPP_1D
796	END SUBROUTINE limit_netcdf
797
798	END MODULE limit
799	!
800	!*******************************************************************************
801

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

Download in other formats:

Original Format