Context Navigation

limit_netcdf.F90 @ 1687

Last change on this file since 1687 was 1687, checked in by Laurent Fairhead, 11 years ago

Nettoyage général pour la bonne marche de makegcm
FH

General housekeeping for makegcm to work as expected
FH

File size: 26.4 KB

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

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

Download in other formats:

Original Format