Context Navigation

datareadnc.F @ 3546

Last change on this file since 3546 was 3382, checked in by tbertrand, 6 months ago
LMDZ.PLUTO Adding more options in newstart from Pluto.old TB
File size: 9.7 KB

Line
1	c=======================================================================
2	SUBROUTINE datareadnc(relief,filename,phisinit,alb,ith,
3	. zmea,zstd,zsig,zgam,zthe)
4	c=======================================================================
5	c
6	c
7	c Author: F. Hourdin 01/1997
8	c -------
9	c
10	c Object: To read data from Martian surface to use in a GCM
11	c ------ from NetCDF file "surface.nc"
12	c
13	c
14	c Arguments:
15	c ----------
16	c
17	c Inputs:
18	c ------
19	c
20	c Outputs:
21	c --------
22	c
23	c=======================================================================
24	c donnees ALBEDO, INERTIE THERMIQUE, RELIEF:
25	c
26	c Ces donnees sont au format NetCDF dans le fichier "surface.nc"
27	c
28	c 360 valeurs en longitude (de -179.5 a 179.5)
29	c 180 valeurs en latitudes (de 89.5 a -89.5)
30	c
31	c Pour les passer au format de la grille, on utilise "interp_horiz.F"
32	c
33	c Il faut donc que ces donnees soient au format grille scalaire
34	c (imold+1 jmold+1)
35	c avec passage des coordonnees de la "boite" (rlonu, rlatv)
36	c
37	c On prend imd (d pour donnees!)
38	c imd = 360 avec copie de la 1ere valeur sur la imd+1
39	c (rlonud de -179 a -181)
40	c jmd = 179
41	c (rlatvd de 89 a -89)
42	c=======================================================================
43
44	use datafile_mod, only: datadir, surfdir
45	! to use 'getin'
46	USE ioipsl_getincom
47	USE comconst_mod, ONLY: g,pi
48	implicit none
49
50	#include "dimensions.h"
51	#include "paramet.h"
52	#include "comgeom.h"
53	#include "netcdf.inc"
54
55	c=======================================================================
56	c Declarations:
57	C=======================================================================
58
59	INTEGER imd,jmd,imdp1,jmdp1
60	parameter (imd=1440,jmd=719,imdp1=1441,jmdp1=720)
61
62	INTEGER iimp1
63	parameter (iimp1=iim+1-1/iim)
64
65	character(len=3),intent(inout) :: relief*3
66	character(len=*),intent(in) :: filename ! surface.nc file
67	real,intent(out) :: phisinit(iimp1*jjp1) ! surface geopotential
68	real,intent(out) :: alb(iimp1*jjp1) ! albedo
69	real,intent(out) :: ith(iimp1*jjp1) ! thermal inertia
70	real,intent(out) :: zmea(imdp1*jmdp1)
71	real,intent(out) :: zstd(imdp1*jmdp1)
72	real,intent(out) :: zsig(imdp1*jmdp1)
73	real,intent(out) :: zgam(imdp1*jmdp1)
74	real,intent(out) :: zthe(imdp1*jmdp1)
75
76	REAL zdata(imd*jmdp1)
77	REAL zdataS(imdp1*jmdp1)
78	REAL pfield(iimp1*jjp1)
79
80	INTEGER ierr
81
82	INTEGER unit,nvarid
83
84	INTEGER i,j,k
85
86	INTEGER klatdat,ngridmixgdat
87	PARAMETER (klatdat=720,ngridmixgdat=1440)
88
89	c on passe une grille en rlonu rlatv et im+1 jm a interp_horiz)
90
91	REAL longitude(imd),latitude(jmdp1) ! Pour lecture des donnees
92	REAL rlonud(imdp1),rlatvd(jmd)
93
94	CHARACTER*20 string
95	DIMENSION string(4)
96	!#include "fxyprim.h"
97
98	pi=2.*ASIN(1.)
99
100	c=======================================================================
101	c rlonud, rlatvd
102	c=======================================================================
103
104	c-----------------------------------------------------------------------
105	c Lecture NetCDF des donnees latitude et longitude
106	c-----------------------------------------------------------------------
107	ierr = NF_OPEN (trim(datadir)//'/'//trim(surfdir)//'/'//
108	& trim(adjustl(filename)),
109	& NF_NOWRITE,unit)
110	IF (ierr.NE.NF_NOERR) THEN
111	! In ye old days this file was stored in datadir;
112	! let's be retro-compatible
113	ierr = NF_OPEN (trim(datadir)//'/'//
114	& trim(adjustl(filename)),
115	& NF_NOWRITE,unit)
116
117	ENDIF
118	IF (ierr.NE.NF_NOERR) THEN
119	write(,)'Error : cannot open file '//trim(filename)
120	write(,)'(in phystd/datareadnc.F)'
121	write(,)'It should be in :',trim(datadir),'/',trim(surfdir)
122	write(,)'Check that your path to datagcm:',trim(datadir)
123	write(,)' is correct. You can change it in callphys.def with:'
124	write(,)' datadir = /absolute/path/to/datagcm'
125	write(,)'If necessary surface.nc (and other datafiles)'
126	write(,)' can be obtained online on:'
127	write(,)' http://www.lmd.jussieu.fr/~lmdz/planets/'//
128	& 'LMDZ.GENERIC/datagcm/'
129	STOP
130	ENDIF
131
132	c
133	c Lecture des latitudes (coordonnees):
134	c
135	ierr = NF_INQ_VARID (unit, "latitude", nvarid)
136	#ifdef NC_DOUBLE
137	ierr = NF_GET_VAR_DOUBLE(unit, nvarid, latitude)
138	#else
139	ierr = NF_GET_VAR_REAL(unit, nvarid, latitude)
140	#endif
141	c
142	c Lecture des longitudes (coordonnees):
143	c
144	ierr = NF_INQ_VARID (unit, "longitude", nvarid)
145	#ifdef NC_DOUBLE
146	ierr = NF_GET_VAR_DOUBLE(unit, nvarid, longitude)
147	#else
148	ierr = NF_GET_VAR_REAL(unit, nvarid, longitude)
149	#endif
150
151	c-----------------------------------------------------------------------
152	c Passage au format boites scalaires
153	c-----------------------------------------------------------------------
154
155	c-----------------------------------------------------------------------
156	c longitude(imd) --> rlonud(imdp1)
157	c-----------------------------------------------------------------------
158
159	c Passage en coordonnees boites scalaires et en radian
160	do i=1,imd
161	rlonud(i)=(longitude(i)+.5)*pi/180.
162	enddo
163
164	c Repetition de la valeur im+1
165	rlonud(imdp1)=rlonud(1) + 2*pi
166
167	c-----------------------------------------------------------------------
168	c latitude(jmdp1) --> rlonvd(jmd)
169	c-----------------------------------------------------------------------
170
171	c Passage en coordonnees boites scalaires et en radian
172	do j=1,jmd
173	rlatvd(j)=(latitude(j)-.5)*pi/180.
174	enddo
175
176	c=======================================================================
177	c lecture NetCDF de albedo, thermal, relief, zdtm (pour francois Lott)
178	c=======================================================================
179
180	string(1) = 'albedo'
181	string(2) = 'thermal'
182	if (relief.ne.'pla') then
183	write(,) ' La topographie est celle de MOLA'
184	relief = 'MOL'
185	string(3) = 'z'//relief
186	else
187	string(3) = 'zMOL' ! pour qu''il lise qqchose sur le fichier
188	! remise a 0 derriere
189	endif
190	string(4) = 'zMOL' ! lecture pour calcul topog. sous-maille
191
192
193	DO k=1,4
194	write(,) 'string',k,string(k)
195
196	c-----------------------------------------------------------------------
197	c initialisation
198	c-----------------------------------------------------------------------
199	pfield(1:iimp1*jjp1)=0
200	zdata(1:imd*jmdp1)=0
201	zdataS(1:iimp1*jjp1)=0
202
203	c-----------------------------------------------------------------------
204	c Lecture NetCDF
205	c-----------------------------------------------------------------------
206
207	ierr = NF_INQ_VARID (unit, string(k), nvarid)
208	#ifdef NC_DOUBLE
209	ierr = NF_GET_VAR_DOUBLE(unit, nvarid, zdata)
210	#else
211	ierr = NF_GET_VAR_REAL(unit, nvarid, zdata)
212	#endif
213
214	c-----------------------------------------------------------------------
215	c Cas particulier "Francois Lott" ( k=4 ) (relief sous-maille)
216	c-----------------------------------------------------------------------
217	if (k.eq.4) then
218
219	zdata(:)=1000.*zdata(:)
220	longitude(:)=(pi/180.)*longitude(:)
221	latitude(:)=(pi/180.)*latitude(:)
222
223	call grid_noro1(1440,720, longitude, latitude, zdata,
224	. iim, jjp1, rlonv, rlatu, zmea,zstd,zsig,zgam,zthe)
225
226	!CALL dump2d(iip1,jjp1,zmea,'zmea')
227	!CALL dump2d(iip1,jjp1,zstd,'zstd')
228	!CALL dump2d(iip1,jjp1,zsig,'zsig')
229	!CALL dump2d(iip1,jjp1,zgam,'zgam')
230	!CALL dump2d(iip1,jjp1,zthe,'zthe')
231
232	endif
233
234	c-----------------------------------------------------------------------
235	c Passage de zdata en grille (imdp1 jmdp1)
236	c-----------------------------------------------------------------------
237	do j=1,jmdp1
238	do i=1,imd
239	zdataS(i+imdp1(j-1)) = zdata(i+ngridmixgdat(j-1))
240	enddo
241	zdataS(imdp1+imdp1(j-1)) = zdata(1+ngridmixgdat(j-1))
242	enddo
243
244	c-----------------------------------------------------------------------
245	c Interpolation
246	c-----------------------------------------------------------------------
247	call interp_horiz(zdataS,pfield,imd,jmd,
248	. iim, jjm,1,rlonud,rlatvd,rlonu,rlatv)
249
250	c-----------------------------------------------------------------------
251	c Periodicite
252	c-----------------------------------------------------------------------
253
254	do j=1,jjp1
255	pfield(iimp1j) = pfield(1+iimp1(j-1))
256	enddo
257
258	c-----------------------------------------------------------------------
259	c Sauvegarde des champs
260	c-----------------------------------------------------------------------
261
262	if (k.eq.1) then ! albedo
263	do i=1,iimp1*jjp1
264	alb(i) = pfield(i)
265	enddo
266	elseif (k.eq.2) then ! thermal
267	do i=1,iimp1*jjp1
268	ith(i) = pfield(i)
269	enddo
270	elseif (k.eq.3) then ! relief
271	if (relief.eq.'pla') then
272	phisinit(1:iimp1*jjp1)=0
273	else
274	phisinit(1:iimp1jjp1)=pfield(1:iimp1jjp1)
275	endif
276	endif
277
278	ENDDO
279
280	c-----------------------------------------------------------------------
281	c Traitement Phisinit
282	c-----------------------------------------------------------------------
283
284	phisinit(1:iimp1jjp1)=1000.phisinit(1:iimp1*jjp1)
285	!CALL dump2d(iimp1,jjp1,phisinit,'Altitude en m')
286	phisinit(:)=g*phisinit(:)
287
288	c-----------------------------------------------------------------------
289	c FIN
290	c-----------------------------------------------------------------------
291
292	END

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

Download in other formats:

Original Format