Context Navigation

datareadnc.F @ 601

Last change on this file since 601 was 588, checked in by emillour, 13 years ago

Generic GCM:
Some cleanup and bug fixing:

"cloudfrac" was not well written to restartfi (wrong size).
missing save attribute for "reffrad" in physiq.F90.
cleanup recomputation of surface pressure in newstart and change loop order in interp_horiz (which "fixes" an odd behaviour which fills some arrays with zeros, but only when using some versions of ifort!)

File size: 9.7 KB

Line
1	c=======================================================================
2	SUBROUTINE datareadnc(relief,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
45	! to use 'getin'
46	USE ioipsl_getincom
47	implicit none
48
49	#include "dimensions.h"
50	#include "paramet.h"
51	#include "comgeom.h"
52	#include "comconst.h"
53	#include "netcdf.inc"
54
55	c=======================================================================
56	c Declarations:
57	C=======================================================================
58
59	INTEGER imd,jmd,imdp1,jmdp1
60	parameter (imd=360,jmd=179,imdp1=361,jmdp1=180)
61
62	INTEGER iimp1
63	parameter (iimp1=iim+1-1/iim)
64
65	CHARACTER relief*3
66
67	REAL zdata(imd*jmdp1)
68	REAL zdataS(imdp1*jmdp1)
69	REAL pfield(iimp1*jjp1)
70
71	REAL alb(iimp1*jjp1)
72	REAL ith(iimp1*jjp1)
73	REAL phisinit(iimp1*jjp1)
74
75	REAL zmea(imdp1*jmdp1)
76	REAL zstd(imdp1*jmdp1)
77	REAL zsig(imdp1*jmdp1)
78	REAL zgam(imdp1*jmdp1)
79	REAL zthe(imdp1*jmdp1)
80
81	INTEGER ierr
82
83	INTEGER unit,nvarid
84
85	INTEGER i,j,k
86
87	INTEGER klatdat,ngridmxgdat
88	PARAMETER (klatdat=180,ngridmxgdat=360)
89
90	c on passe une grille en rlonu rlatv et im+1 jm a interp_horiz)
91
92	REAL longitude(imd),latitude(jmdp1) ! Pour lecture des donnees
93	REAL rlonud(imdp1),rlatvd(jmd)
94
95	CHARACTER*20 string
96	DIMENSION string(4)
97
98	CHARACTER*20 filename
99	CHARACTER*50 filestring
100
101	#include "lmdstd.h"
102	#include "fxyprim.h"
103
104	pi=2.*ASIN(1.)
105
106	c=======================================================================
107	c rlonud, rlatvd
108	c=======================================================================
109
110	c-----------------------------------------------------------------------
111	c Lecture NetCDF des donnees latitude et longitude
112	c-----------------------------------------------------------------------
113	! First get the corret value of dtadir, if not already done:
114	! default 'datadir' is set in "datadir_mod"
115	call getin("datadir",datadir)
116	write(,) 'Choice of surface data is:'
117	filestring='ls '//trim(datadir)//'/*.nc'
118	call system(filestring)
119
120	write(,) ''
121	write(,) 'Please choose the relevant file:'
122	write(,) '(e.g. type "surface_mars.nc")'
123	read(*,fmt='(a20)') filename
124
125	ierr = NF_OPEN (trim(datadir)//'/'//trim(adjustl(filename)),
126	& NF_NOWRITE,unit)
127	IF (ierr.NE.NF_NOERR) THEN
128	write(,)'Error : cannot open file '//trim(filename)
129	write(,)'(in phystd/datareadnc.F)'
130	write(,)'It should be in :',trim(datadir),'/'
131	write(,)'Check that your path to datagcm:',trim(datadir)
132	write(,)' is correct. You can change it in callphys.def with:'
133	write(,)' datadir = /absolute/path/to/datagcm'
134	write(,)'If necessary surface.nc (and other datafiles)'
135	write(,)' can be obtained online on:'
136	write(,)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
137	CALL ABORT
138	ENDIF
139
140	c
141	c Lecture des latitudes (coordonnees):
142	c
143	ierr = NF_INQ_VARID (unit, "latitude", nvarid)
144	#ifdef NC_DOUBLE
145	ierr = NF_GET_VAR_DOUBLE(unit, nvarid, latitude)
146	#else
147	ierr = NF_GET_VAR_REAL(unit, nvarid, latitude)
148	#endif
149	c
150	c Lecture des longitudes (coordonnees):
151	c
152	ierr = NF_INQ_VARID (unit, "longitude", nvarid)
153	#ifdef NC_DOUBLE
154	ierr = NF_GET_VAR_DOUBLE(unit, nvarid, longitude)
155	#else
156	ierr = NF_GET_VAR_REAL(unit, nvarid, longitude)
157	#endif
158
159	c-----------------------------------------------------------------------
160	c Passage au format boites scalaires
161	c-----------------------------------------------------------------------
162
163	c-----------------------------------------------------------------------
164	c longitude(imd) --> rlonud(imdp1)
165	c-----------------------------------------------------------------------
166
167	c Passage en coordonnees boites scalaires et en radian
168	do i=1,imd
169	rlonud(i)=(longitude(i)+.5)*pi/180.
170	enddo
171
172	c Repetition de la valeur im+1
173	rlonud(imdp1)=rlonud(1) + 2*pi
174
175	c-----------------------------------------------------------------------
176	c latitude(jmdp1) --> rlonvd(jmd)
177	c-----------------------------------------------------------------------
178
179	c Passage en coordonnees boites scalaires et en radian
180	do j=1,jmd
181	rlatvd(j)=(latitude(j)-.5)*pi/180.
182	enddo
183
184	c=======================================================================
185	c lecture NetCDF de albedo, thermal, relief, zdtm (pour francois Lott)
186	c=======================================================================
187
188	string(1) = 'albedo'
189	string(2) = 'thermal'
190	if (relief.ne.'pla') then
191	write(,) ' La topographie est celle de MOLA'
192	relief = 'MOL'
193	string(3) = 'z'//relief
194	else
195	string(3) = 'zMOL' ! pour qu''il lise qqchose sur le fichier
196	! remise a 0 derriere
197	endif
198	string(4) = 'zMOL' ! lecture pour calcul topog. sous-maille
199
200
201	DO k=1,4
202	write(,) 'string',k,string(k)
203
204	c-----------------------------------------------------------------------
205	c initialisation
206	c-----------------------------------------------------------------------
207	call initial0(iimp1*jjp1,pfield)
208	call initial0(imd*jmdp1,zdata)
209	call initial0(imdp1*jmdp1,zdataS)
210
211	c-----------------------------------------------------------------------
212	c Lecture NetCDF
213	c-----------------------------------------------------------------------
214
215	ierr = NF_INQ_VARID (unit, string(k), nvarid)
216	#ifdef NC_DOUBLE
217	ierr = NF_GET_VAR_DOUBLE(unit, nvarid, zdata)
218	#else
219	ierr = NF_GET_VAR_REAL(unit, nvarid, zdata)
220	#endif
221
222	c-----------------------------------------------------------------------
223	c Cas particulier "Francois Lott" ( k=4 ) (relief sous-maille)
224	c-----------------------------------------------------------------------
225	if (k.eq.4) then
226
227	call multscal(imd*jmdp1,zdata,1000.,zdata)
228	call multscal(imd,longitude,pi/180.,longitude)
229	call multscal(jmdp1,latitude,pi/180.,latitude)
230
231	call grid_noro1(360, 180, longitude, latitude, zdata,
232	. iim, jjp1, rlonv, rlatu, zmea,zstd,zsig,zgam,zthe)
233
234	CALL dump2d(iip1,jjp1,zmea,'zmea')
235	CALL dump2d(iip1,jjp1,zstd,'zstd')
236	CALL dump2d(iip1,jjp1,zsig,'zsig')
237	CALL dump2d(iip1,jjp1,zgam,'zgam')
238	CALL dump2d(iip1,jjp1,zthe,'zthe')
239
240	endif
241
242	c-----------------------------------------------------------------------
243	c Passage de zdata en grille (imdp1 jmdp1)
244	c-----------------------------------------------------------------------
245	do j=1,jmdp1
246	do i=1,imd
247	zdataS(i+imdp1(j-1)) = zdata(i+ngridmxgdat(j-1))
248	enddo
249	zdataS(imdp1+imdp1(j-1)) = zdata(1+ngridmxgdat(j-1))
250	enddo
251
252	c-----------------------------------------------------------------------
253	c Interpolation
254	c-----------------------------------------------------------------------
255	call interp_horiz(zdataS,pfield,imd,jmd,
256	. iim, jjm,1,rlonud,rlatvd,rlonu,rlatv)
257
258	c-----------------------------------------------------------------------
259	c Periodicite
260	c-----------------------------------------------------------------------
261
262	do j=1,jjp1
263	pfield(iimp1j) = pfield(1+iimp1(j-1))
264	enddo
265
266	c-----------------------------------------------------------------------
267	c Sauvegarde des champs
268	c-----------------------------------------------------------------------
269
270	if (k.eq.1) then ! albedo
271	do i=1,iimp1*jjp1
272	alb(i) = pfield(i)
273	enddo
274	elseif (k.eq.2) then ! thermal
275	do i=1,iimp1*jjp1
276	ith(i) = pfield(i)
277	enddo
278	elseif (k.eq.3) then ! relief
279	if (relief.eq.'pla') then
280	call initial0(iimp1*jjp1,phisinit)
281	else
282	do i=1,iimp1*jjp1
283	phisinit(i) = pfield(i)
284	enddo
285	endif
286	endif
287
288	ENDDO
289
290	c-----------------------------------------------------------------------
291	c Traitement Phisinit
292	c-----------------------------------------------------------------------
293
294	DO i=1,iimp1*jjp1
295	phisinit(i)=1000.*phisinit(i)
296	ENDDO
297	CALL dump2d(iimp1,jjp1,phisinit,'Altitude en m')
298	CALL multscal(iimp1*jjp1,phisinit,g,phisinit)
299
300	c-----------------------------------------------------------------------
301	c FIN
302	c-----------------------------------------------------------------------
303
304	END

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

Download in other formats:

Original Format