Context Navigation

limit_netcdf.F90 @ 2889

Last change on this file since 2889 was 2889, checked in by dcugnet, 7 years ago

Change of behaviour for ce0l:

The time interpolation of AMIP mid-month boundary conditions datasets (corresponding ce0l standard files/variables names: "amipbc_sst_1x1.bnc/tosbcs", "amipbc_sic_1x1.nc/sicbcs") is no longer using splines, but is now linear. See https://pcmdi.llnl.gov/mips/amip/details
New ce0l standard files/variables names "amib_sst_1x1.nc/sst", "amip_sic_1x1.nc/sic" are introduced to allow splines time interpolation with monthly mean observations if needed (not adviced).

Note that splines time interpolation is still used for other fields "BILS", "ALB", "RUG"
and other cases ("amip_sst_1x1.nc" "cpl_atm_sst.nc" "histmth_sst.nc" "sstk.nc" for SST).

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

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

Download in other formats:

Original Format