Context Navigation

globgrd.f90 @ 8379

Last change on this file since 8379 was 8379, checked in by jan.polcher, 5 months ago
A bug which was corrected in the tunk but not 2.2
File size: 40.8 KB

Line
1	! ==============================================================================================================================\n
2	! MODULE globgrd : This module is dedicated to managing the spatial grid of the forcing. It can either read a file
3	! containing the grid information, as is the case for WRF forcing, or obtain the grid from the forcing files.
4	! The module has also the possibility to create a grid description files for certain applications like
5	! for instance in a coupling of ORCHIDEE through OASIS.
6	! For this purpose the module provides 4 subroutines :
7	! globgrd_getdomsz : This routine allows to get the domain size of the forcing based on a file it will explore.
8	! globgrd_getgrid : This routine extracts the coordinates and land/sea mask from the domain files.
9	! globgrd_writevar : Writes a variables into a netCDF file which can then be analysed to verify that the
10	! forcing grid was well read and interpreted by this module.
11	! globgrd_writegrid : Write a grid description file which has the WRF flavor. It allows to exchange grid information
12	! between an atmospheric model (mostly a driver !) and ORCHIDEE which are coupled through OASIS.
13	!
14	! CONTACT : jan.polcher@lmd.jussieu.fr
15	!
16	! LICENCE : IPSL (2016)
17	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
18	!
19	!>\BRIEF
20	!!
21	!! RECENT CHANGE(S): None
22	!!
23	!! REFERENCE(S) : None
24	!!
25	!_ ================================================================================================================================
26	MODULE globgrd
27	!
28	!
29	USE defprec
30	USE netcdf
31	!
32	USE ioipsl
33	!
34	USE grid
35	USE forcing_tools
36	!
37	IMPLICIT NONE
38	!
39	PRIVATE
40	PUBLIC :: globgrd_getdomsz, globgrd_getgrid, globgrd_writevar, globgrd_writegrid
41	!
42	!
43	LOGICAL, SAVE :: is_forcing_file=.FALSE.
44	!
45	CONTAINS
46	!!
47	!! =============================================================================================================================
48	!! SUBROUTINE: globgrd_getdomsz
49	!!
50	!>\BRIEF This routine allows to get the domain size of the forcing based on a file it will explore.
51	!!
52	!! DESCRIPTION: The routine opens the file and explores it. It can either be a forcing file or a grid description
53	!! file from WRF. Progressively this should be opened to other ways of describing the grid over which
54	!! the forcing is provided.
55	!! The routing will return the sizes in I and J and the number of land points.
56	!! The zooming interval is also provided so that only the dimensions over the domain used can be computed.
57	!!
58	!! \n
59	!_ ==============================================================================================================================
60	!!
61	!---------------------------------------------------------------------
62	!-
63	!-
64	!---------------------------------------------------------------------
65	SUBROUTINE globgrd_getdomsz(filename, iim, jjm, nbland, model_guess, fid, forcingfile, zoom_lon, zoom_lat)
66	!
67	! INPUT
68	!
69	CHARACTER(LEN=*), INTENT(in) :: filename
70	CHARACTER(LEN=*), INTENT(in), OPTIONAL :: forcingfile(:)
71	REAL(r_std), DIMENSION(2), INTENT(in), OPTIONAL :: zoom_lon, zoom_lat
72	!
73	! OUTPUT
74	!
75	INTEGER(i_std), INTENT(out) :: fid
76	INTEGER(i_std), INTENT(out) :: iim, jjm, nbland
77	CHARACTER(LEN=*), INTENT(out) :: model_guess
78	!
79	! LOCAL
80	!
81	INTEGER(i_std) :: iret, ndims, nvars, nb_atts, id_unlim, iv, lll
82	INTEGER(i_std) :: iindex_init, jindex_init, iindex_end, jindex_end
83	INTEGER(i_std) :: iim_full, jjm_full, nbland_full
84	CHARACTER(LEN=20) :: axname, varname
85	CHARACTER(LEN=120) :: tmpfile
86	REAL(r_std), DIMENSION(2) :: loczoom_lon, loczoom_lat
87	INTEGER(i_std), DIMENSION(2) :: tmp_lon_ind1, tmp_lat_ind1, tmp_lon_ind2, tmp_lat_ind2
88	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask,zoom_mask, lat, lon
89	!
90	! Set default values against which we can test
91	!
92	iim = -1
93	jjm = -1
94	!
95	! Verify the grid file name
96	! if forcing_file exists and if grid file exists: tmpfile = grid file
97	! if forcing_file exits and grid file not exists: tmpfile = forcing_file
98	! if forcing_file not exists: tmp_file= grid file
99	!
100	IF ( PRESENT(forcingfile) ) THEN
101	is_forcing_file=.TRUE.
102	IF ( INDEX(filename,"NONE") >= 1 ) THEN
103	tmpfile=forcingfile(1)
104	ELSE
105	tmpfile=filename
106	ENDIF
107	ELSE
108	is_forcing_file=.FALSE.
109	tmpfile=filename
110	ENDIF
111	!
112	! Verify that the zoomed region is provided. Else choose the entire globe
113	!
114	IF ( PRESENT(zoom_lon) .AND. PRESENT(zoom_lat) ) THEN
115	loczoom_lon = zoom_lon
116	loczoom_lat = zoom_lat
117	ELSE
118	loczoom_lon(1) = -180.0
119	loczoom_lon(2) = 180.0
120	loczoom_lat(1) = -90.0
121	loczoom_lat(2) = 90.0
122	ENDIF
123	!
124	! Open the correct file
125	!
126	iret = NF90_OPEN (tmpfile, NF90_NOWRITE, fid)
127	IF (iret /= NF90_NOERR) THEN
128	CALL ipslerr (3,'globgrd_getdomsz',"Error opening the grid file :",tmpfile, " ")
129	ENDIF
130	!
131	!
132	iret = NF90_INQUIRE (fid, nDimensions=ndims, nVariables=nvars, &
133	nAttributes=nb_atts, unlimitedDimId=id_unlim)
134	IF (iret /= NF90_NOERR) THEN
135	CALL ipslerr (3,'globgrd_getdomsz',"Error in NF90_INQUIRE :",tmpfile, " ")
136	ENDIF
137	!
138	DO iv=1,ndims
139	!
140	iret = NF90_INQUIRE_DIMENSION (fid, iv, name=axname, len=lll)
141	IF (iret /= NF90_NOERR) THEN
142	CALL ipslerr (3,'globgrd_getdomsz',"Could not get size of dimension :"," "," ")
143	ENDIF
144	!
145	! This can be refined by testing the actual grid found in the file.
146	!
147	SELECT CASE(axname)
148	!
149	!! Coordinate variables used by WRF.
150	CASE("west_east")
151	iim_full = lll
152	model_guess = "WRF"
153	CASE("south_north")
154	jjm_full = lll
155	model_guess = "WRF"
156	!
157	!! Variables used in WFDEI
158	CASE("lon")
159	iim_full = lll
160	model_guess = "regular"
161	CASE("lat")
162	jjm_full = lll
163	model_guess = "regular"
164	CASE("nbland")
165	nbland_full = lll
166	!
167	!! Variables used by CRU-NCEP
168	CASE("nav_lon")
169	iim_full = lll
170	model_guess = "regular"
171	CASE("nav_lat")
172	jjm_full = lll
173	model_guess = "regular"
174	CASE("land")
175	nbland_full = lll
176
177
178	!! Variables used by CRU-JRA (v2.2.2)
179	CASE("longitude")
180	iim_full = lll
181	model_guess = "regular"
182	CASE("latitude")
183	jjm_full = lll
184	model_guess = "regular"
185	END SELECT
186	ENDDO
187	!
188	! If we have a WRF file we need to count the number of land points, define iim jjm and mask
189	!
190	IF ( model_guess == "WRF" ) THEN
191
192	IF ( .NOT. ALLOCATED(mask) ) ALLOCATE(mask(iim_full,jjm_full))
193
194	varname = "LANDMASK"
195	iret = NF90_INQ_VARID (fid, varname, iv)
196	IF (iret /= NF90_NOERR) THEN
197	CALL ipslerr (3,'globgrd_getdomsz',"Could not find variable ", varname," ")
198	ELSE
199	iret = NF90_GET_VAR (fid,iv,mask)
200	ENDIF
201
202	nbland_full = COUNT(mask > 0.5)
203
204
205	IF ( .NOT. ALLOCATED(lat)) ALLOCATE(lat(iim_full,jjm_full))
206	IF ( .NOT. ALLOCATED(lon)) ALLOCATE(lon(iim_full,jjm_full))
207
208	varname = "XLONG_M"
209	iret = NF90_INQ_VARID (fid, varname, iv)
210	IF (iret /= NF90_NOERR) THEN
211	CALL ipslerr (3,'globgrd_getdomsz',"Could not find variable ", varname," ")
212	ELSE
213	iret = NF90_GET_VAR (fid,iv,lon)
214	ENDIF
215
216	varname = "XLAT_M"
217	iret = NF90_INQ_VARID (fid, varname, iv)
218	IF (iret /= NF90_NOERR) THEN
219	CALL ipslerr (3,'globgrd_getdomsz',"Could not find variable ", varname," ")
220	ELSE
221	iret = NF90_GET_VAR (fid,iv,lat)
222	ENDIF
223
224	!define nbland for full-region simulation in case of need
225	nbland = nbland_full
226
227	ENDIF
228	!
229	!
230	! If we are in the case of a forcing file and model_guess is regular grid,
231	! then a few functions from forcing_tools need to be called
232	! so that the file is analysed with the tools of the forcing module.
233	!
234	! predefine nbland, iim, jjm in the case is_forcing_file=false, or full grid simulation for wrf
235	iim = iim_full
236	jjm = jjm_full
237	IF ( is_forcing_file .AND. model_guess .EQ. "regular") THEN
238	!
239	! Because we are re-using routines from the forcing module, we have to
240	! close the file. It will be opened again by the forcing module.
241	!
242	iret = NF90_CLOSE(fid)
243	IF (iret /= NF90_NOERR) THEN
244	CALL ipslerr (3,'globgrd_getdomzz',"Error closing the output file :",filename, " ")
245	ENDIF
246	!
247	! Set last argument closefile=.FALSE. as the forcing file has been closed here above.
248	! This will also induce that dump_mask=.FALSE. in forcing_getglogrid and the
249	! file forcing_mask_glo.nc will not be created. See also ticket #691
250	CALL forcing_getglogrid (1, forcingfile, iim_full, jjm_full, nbland_full, .FALSE.)
251	WRITE(,) forcingfile, "Forcing file with dimensions : ", iim_full, jjm_full, nbland_full
252	!
253	CALL forcing_zoomgrid (loczoom_lon, loczoom_lat, forcingfile(1), model_guess, .TRUE.)
254	!
255	CALL forcing_givegridsize (iim, jjm, nbland)
256	!
257
258	ELSE IF ( is_forcing_file .AND. model_guess .EQ. "WRF") THEN
259	!
260	! if forcing_file exists and model_guess is wrf, and zoomed domain is asked,
261	! then we define the domain size according to the zoom
262	! if the full domain is asked, then we do not define the domain size again here.
263	!
264	! get the beginning and endding index in the case of zoomed region, for WRF grids (XW)
265	! Not to close the grid file here, to be used by globgrd_getgrid
266	!
267	!IF ( PRESENT(zoom_lon) .AND. PRESENT(zoom_lat) ) THEN
268	IF ( PRESENT(zoom_lon) .OR. PRESENT(zoom_lat) ) THEN
269
270	! to get new iim, jjm and nbland for the zoomed region
271	tmp_lon_ind1 = MINLOC(ABS(lon(:,:)-loczoom_lon(1)))
272	tmp_lat_ind1 = MINLOC(ABS(lat(:,:)-loczoom_lat(1)))
273	tmp_lon_ind2 = MINLOC(ABS(lon(:,:)-loczoom_lon(2)))
274	tmp_lat_ind2 = MINLOC(ABS(lat(:,:)-loczoom_lat(2)))
275	!
276	! the zoomed region
277	! lon1, lat2------- lon2,lat2
278	! \| \|
279	! lon1, lat1------- lon2,lat1
280	!
281	iindex_init = tmp_lon_ind1(1)
282	jindex_init= tmp_lat_ind1(2)
283
284	iindex_end= tmp_lon_ind2(1)
285	jindex_end= tmp_lat_ind2(2)
286
287	iim = iindex_end - iindex_init + 1
288	jjm = jindex_end - jindex_init + 1
289
290
291	IF ( .NOT. ALLOCATED(zoom_mask) ) ALLOCATE(zoom_mask(iim,jjm))
292	zoom_mask = mask(iindex_init:iindex_end, jindex_init:jindex_end)
293	nbland = COUNT(zoom_mask > 0.5)
294	IF ( ALLOCATED(zoom_mask) ) DEALLOCATE(zoom_mask)
295
296	ENDIF
297
298	IF ( ALLOCATED(lat) ) DEALLOCATE(lat)
299	IF ( ALLOCATED(lon) ) DEALLOCATE(lon)
300	IF ( ALLOCATED(mask) ) DEALLOCATE(mask)
301
302	ENDIF
303	!
304	! Do a final test to see if we got the information needed.
305	!
306	IF ( iim < 0 .OR. jjm < 0 ) THEN
307	CALL ipslerr (3,'globgrd_getdomsz',"Could not get the horizontal size of the domaine out of the file",&
308	& filename,"Are you sure that the case for this type of file is foreseen ? ")
309	ENDIF
310	!
311	!
312	END SUBROUTINE globgrd_getdomsz
313	!!
314	!! =============================================================================================================================
315	!! SUBROUTINE: globgrd_getgrid
316	!!
317	!>\BRIEF This routine extracts the coordinates and land/sea mask from the domain files.
318	!!
319	!! DESCRIPTION: The domain size is provided together with the netCDF file ID so that the main information can
320	!! be extracted from the file. We will read the longitude, latitude, land/sea mask and calendar.
321	!! This allows to set-up ORCHIDEE. We also provide the corners of the grid-boxes as this is needed
322	!! for setting-up OASIS but is computed more correctly in grid.f90 for ORCHIDEE.
323	!! This routine is only an interface to globgrd_getwrf, globgrd_getregular and forcing_givegrid.
324	!! forcing_givegrid is an interface to the forcing_tools.f90 module so that we are certain to have
325	!! the same grid information between both modules.
326	!!
327	!! \n
328	!_ ==============================================================================================================================
329	!!
330	!---------------------------------------------------------------------
331	!-
332	!-
333	!---------------------------------------------------------------------
334	SUBROUTINE globgrd_getgrid(fid, iim, jjm, nbland, model_guess, lon, lat, mask, area, corners, &
335	& lindex, contfrac, calendar, zoom_lon, zoom_lat)
336	!
337	!
338	! This subroutine only switched between routines to extract and compte the grid data needed for
339	! ORCHIDEE.
340	!
341	!
342	! INPUT
343	!
344	INTEGER(i_std), INTENT(in) :: fid
345	INTEGER(i_std), INTENT(in) :: iim, jjm, nbland
346	CHARACTER(LEN=*), INTENT(in) :: model_guess
347	REAL(r_std), DIMENSION(2), INTENT(in), OPTIONAL :: zoom_lon, zoom_lat
348	!
349	! OUTPUT
350	!
351	REAL(r_std),DIMENSION(iim,jjm), INTENT(out) :: lon, lat, mask, area
352	REAL(r_std),DIMENSION(iim,jjm,4,2), INTENT(out) :: corners
353	INTEGER(i_std), DIMENSION(nbland), INTENT(out) :: lindex
354	REAL(r_std),DIMENSION(nbland), INTENT(out) :: contfrac
355	CHARACTER(LEN=20), INTENT(out) :: calendar
356	!
357	SELECT CASE(model_guess)
358
359	CASE("WRF")
360	CALL globgrd_getwrf(fid, iim, jjm, nbland, lon, lat, mask, area, corners, &
361	& lindex, contfrac, calendar, zoom_lon, zoom_lat)
362	CASE("regular")
363	IF ( .NOT. is_forcing_file ) THEN
364	CALL globgrd_getregular(fid, iim, jjm, nbland, lon, lat, mask, area, corners, &
365	& lindex, contfrac, calendar)
366	ELSE
367	CALL forcing_givegrid(lon, lat, mask, area, corners, lindex, contfrac, calendar)
368	CALL forcing_close()
369	ENDIF
370	CASE DEFAULT
371	CALL ipslerr (3,'globgrd_getgrid',"The model/grid type we guessed is not recognized here. model_guess =",&
372	& model_guess,"Have you used the right file and are you sure that this case is foreseen ? ")
373	END SELECT
374	!
375	!
376	END SUBROUTINE globgrd_getgrid
377	!!
378	!! =============================================================================================================================
379	!! SUBROUTINE: globgrd_regular
380	!!
381	!>\BRIEF The routine to obtain regular grids from the file.
382	!!
383	!! DESCRIPTION: Read the regular grid and its information from the opened file (fid).
384	!!
385	!! \n
386	!_ ==============================================================================================================================
387	!!
388	SUBROUTINE globgrd_getregular(fid, iim, jjm, nbland, lon, lat, mask, area, corners, &
389	& lindex, contfrac, calendar)
390	!
391	USE defprec
392	USE netcdf
393	!
394	! INPUT
395	!
396	INTEGER(i_std), INTENT(in) :: fid
397	INTEGER(i_std), INTENT(in) :: iim, jjm, nbland
398	!
399	! OUTPUT
400	!
401	REAL(r_std),DIMENSION(iim,jjm), INTENT(out) :: lon, lat, mask, area
402	REAL(r_std),DIMENSION(iim,jjm,4,2), INTENT(out) :: corners
403	INTEGER(i_std), DIMENSION(nbland), INTENT(out) :: lindex
404	REAL(r_std),DIMENSION(nbland), INTENT(out) :: contfrac
405	CHARACTER(LEN=20), INTENT(out) :: calendar
406	!
407	! LOCAL
408	!
409	INTEGER(i_std) :: iret, iv, nvars, varndim
410	INTEGER(i_std) :: i, j
411	CHARACTER(LEN=20) :: varname
412	INTEGER(i_std), DIMENSION(4) :: vardims
413	REAL(r_std) :: dx
414	!
415	! Set some default values agains which we can check
416	!
417	lon(:,:) = val_exp
418	lat(:,:) = val_exp
419	mask(:,:) = val_exp
420	area(:,:) = val_exp
421	corners(:,:,:,:) = val_exp
422	!
423	lindex(:) = INT(val_exp)
424	contfrac(:) = val_exp
425	!
426	! Get the global attributes from grid file
427	!
428	iret = NF90_GET_ATT(fid, NF90_GLOBAL, 'calendar', calendar)
429	IF (iret /= NF90_NOERR) THEN
430	CALL ipslerr (3,'globgrd_getregular',"Could not read the calendar in grid file.", " ", " ")
431	ENDIF
432	!
433	iret = NF90_INQUIRE (fid, nVariables=nvars)
434	!
435	DO iv = 1,nvars
436	!
437	iret = NF90_INQUIRE_VARIABLE(fid, iv, name=varname, ndims=varndim, dimids=vardims)
438	!
439	!
440	SELECT CASE(varname)
441	!
442	CASE("longitude")
443	IF (varndim == 1 ) THEN
444	DO j=1,jjm
445	iret = NF90_GET_VAR(fid, iv, lon(:,j))
446	ENDDO
447	ELSE IF (varndim == 2 ) THEN
448	iret = NF90_GET_VAR(fid, iv, lon)
449	ELSE
450	CALL ipslerr (3,'globgrd_getregular',"Longitude cannot have more than 2 dimensions","","")
451	ENDIF
452
453	CASE ("latitude")
454	IF (varndim == 1 ) THEN
455	DO i=1,iim
456	iret = NF90_GET_VAR(fid, iv, lat(i,:))
457	ENDDO
458	ELSE IF (varndim == 2 ) THEN
459	iret = NF90_GET_VAR(fid, iv, lon)
460	ELSE
461	CALL ipslerr (3,'globgrd_getregular',"Latitude cannot have more than 2 dimensions","","")
462	ENDIF
463
464	CASE ("mask")
465	IF (varndim /= 2 ) THEN
466	CALL ipslerr (3,'globgrd_getregular',"mask needs to have 2 dimensions","","")
467	ELSE
468	iret = NF90_GET_VAR (fid, iv, mask)
469	ENDIF
470
471	CASE ("areas")
472	IF (varndim /= 2 ) THEN
473	CALL ipslerr (3,'globgrd_getregular',"Areas needs to have 2 dimensions","","")
474	ELSE
475	iret = NF90_GET_VAR (fid, iv, area)
476	ENDIF
477
478	CASE ("corners")
479	IF (varndim /= 4 ) THEN
480	CALL ipslerr (3,'globgrd_getregular',"corners needs to have 4 dimensions","","")
481	ELSE
482	iret = NF90_GET_VAR (fid, iv, corners)
483	ENDIF
484
485	CASE ("landindex")
486	IF (varndim /= 1 ) THEN
487	CALL ipslerr (3,'globgrd_getregular',"landindex is the list of continental points to be gathered", &
488	& "Thus it can only have 1 dimensions","")
489	ELSE
490	iret = NF90_GET_VAR (fid, iv, lindex)
491	ENDIF
492
493	CASE ("contfrac")
494	IF (varndim /= 1 ) THEN
495	CALL ipslerr (3,'globgrd_getregular',"Contfrac needs to be a gathered variable", &
496	& "thus it needs only 1 dimensions","")
497	ELSE
498	iret = NF90_GET_VAR (fid, iv, contfrac)
499	ENDIF
500
501	END SELECT
502	!
503	ENDDO
504	!
505	!
506	iret = NF90_CLOSE(fid)
507	!
508	! Verify that we have al the variables needed to describe the ORCHIDEE grid
509	!
510	IF ( ANY( lon(:,:) == val_exp ) ) THEN
511	CALL ipslerr (3,'globgrd_getregular',"The longitude of the ORCHIDEE grid could not be extracted from the",&
512	& "grid definition file","")
513	ENDIF
514	!
515	IF ( ANY( lat(:,:) == val_exp ) ) THEN
516	CALL ipslerr (3,'globgrd_getregular',"The latitude of the ORCHIDEE grid could not be extracted from the",&
517	& "grid definition file","")
518	ENDIF
519	!
520	IF ( ANY( lindex(:) == INT(val_exp) ) ) THEN
521	CALL ipslerr (3,'globgrd_getregular',"The lindex of the ORCHIDEE grid could not be extracted from the",&
522	& "grid definition file","")
523	ENDIF
524	!
525	IF ( ALL( mask(:,:) == val_exp ) ) THEN
526	CALL ipslerr (3,'globgrd_getregular',"The land mask of the ORCHIDEE grid could not be extracted from the",&
527	& "grid definition file","")
528	ELSE IF (MAXVAL(mask) > 1 ) THEN
529	CALL ipslerr (2,'globgrd_getregular',"We have a special case for the mask which needs to be treated.",&
530	& "The field contains the indices of the land points on a compressed grid.","So we replace them with 1 or 0.")
531	DO i=1,iim
532	DO j=1,jjm
533	IF ( mask(i,j) > iim*jjm ) THEN
534	mask(i,j) = 0
535	ELSE
536	mask(i,j) = 1
537	ENDIF
538	ENDDO
539	ENDDO
540	ENDIF
541	!
542	IF ( ANY( contfrac(:) == val_exp ) ) THEN
543	CALL ipslerr (2,'globgrd_getregular',"The continental fraction of the ORCHIDEE grid could not be extracted from the",&
544	& "grid definition file","Thus on all land points it is set to 1.")
545	contfrac(:) = 1.
546	ENDIF
547	!
548	IF ( ANY( corners(:,:,:,:) == val_exp ) ) THEN
549	CALL ipslerr (3,'globgrd_getregular',"The corners for the ORCHIDEE grid could not be extracted from the",&
550	& "grid definition file","As we have to assume a very general grid we cannot do anything !")
551	ENDIF
552	!
553	!
554	END SUBROUTINE globgrd_getregular
555	!!
556	!! =============================================================================================================================
557	!! SUBROUTINE: globgrd_getwrf
558	!!
559	!>\BRIEF Routine to read the WRF grid description file.
560	!!
561	!! DESCRIPTION: Read the WRF grid and its information from the opened file (fid) and convert
562	!! it to the variables needed by ORCHIDEE.
563	!!
564	!! \n
565	!_ ==============================================================================================================================
566	!!
567	SUBROUTINE globgrd_getwrf(fid, iim, jjm, nbland, lon, lat, mask, area, corners, &
568	& lindex, contfrac, calendar, zoom_lon, zoom_lat)
569	!
570	USE defprec
571	USE netcdf
572	!
573	! INPUT
574	!
575	INTEGER(i_std), INTENT(in) :: fid
576	INTEGER(i_std), INTENT(in) :: iim, jjm, nbland
577	REAL(r_std), DIMENSION(2), INTENT(in), OPTIONAL :: zoom_lon, zoom_lat
578	!
579	! OUTPUT
580	!
581	REAL(r_std),DIMENSION(iim,jjm), INTENT(out) :: lon, lat, mask, area
582	REAL(r_std),DIMENSION(iim,jjm,4,2), INTENT(out) :: corners
583	INTEGER(i_std), DIMENSION(nbland), INTENT(out) :: lindex
584	REAL(r_std),DIMENSION(nbland), INTENT(out) :: contfrac
585	CHARACTER(LEN=20), INTENT(out) :: calendar
586	!
587	! LOCAL
588	!
589	INTEGER(i_std) :: i, ip, jp, k, iret, iv, nvars, varndim
590	INTEGER(i_std),DIMENSION(2) :: tmp_lon_ind1, tmp_lat_ind1, tmp_lon_ind2, tmp_lat_ind2
591	REAL(r_std), DIMENSION(2) :: loczoom_lon, loczoom_lat
592	INTEGER(i_std) :: ndims, nb_atts, id_unlim, lll, iim_full, jjm_full, iim_tmp, jjm_tmp
593	INTEGER(i_std) :: iindex_init, jindex_init, iindex_end, jindex_end, i_orig, j_orig
594	CHARACTER(LEN=20) :: varname
595
596	INTEGER(i_std), DIMENSION(4) :: vardims
597	INTEGER(i_std), DIMENSION(8) :: rose
598	!
599	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask_full, lon_full, lat_full, area_full
600	REAL(r_std),ALLOCATABLE, DIMENSION(:,:) :: mapfac_x_full, mapfac_y_full
601	CHARACTER(LEN=20) :: axname
602	REAL(r_std) :: dx, dy, coslat
603	!
604	REAL(r_std), PARAMETER :: mincos = 0.0001
605	REAL(r_std), PARAMETER :: pi = 3.141592653589793238
606	REAL(r_std), PARAMETER :: R_Earth = 6378000.
607	!
608	! A lot of modifications are added to this subroutine (XW)
609	!
610	! Set some default values agains which we can check
611	!
612	lon(:,:) = val_exp
613	lat(:,:) = val_exp
614	mask(:,:) = val_exp
615	area(:,:) = val_exp
616	corners(:,:,:,:) = val_exp
617	!
618	lindex(:) = INT(val_exp)
619	contfrac(:) = val_exp
620	!
621	calendar = 'gregorian'
622	!
623	! get dimension of full-grid data
624	!
625	iret = NF90_INQUIRE (fid, nDimensions=ndims, nVariables=nvars, &
626	nAttributes=nb_atts, unlimitedDimId=id_unlim)
627	!
628	DO iv=1,ndims
629	!
630	iret = NF90_INQUIRE_DIMENSION (fid, iv, name=axname, len=lll)
631	IF (iret /= NF90_NOERR) THEN
632	CALL ipslerr (3,'globgrd_getdomsz',"Could not get size of dimension :"," "," ")
633	ENDIF
634	!
635	! This can be refined by testing the actual grid found in the file.
636	!
637	SELECT CASE(axname)
638	!
639	!! Coordinate variables used by WRF.
640	CASE("west_east")
641	iim_full = lll
642
643	CASE("south_north")
644	jjm_full = lll
645
646	END SELECT
647	ENDDO
648
649	!
650	! Init projection in grid.f90 so that it can be used later for projections.
651	!
652	CALL grid_initproj(fid, iim_full, jjm_full)
653	!
654	iret = NF90_INQUIRE (fid, nVariables=nvars)
655	!
656	IF (iret /= NF90_NOERR) THEN
657	CALL ipslerr (3,'globgrd_getwrf',"Error inquiering variables from WRF grid file."," ", " ")
658	ENDIF
659	IF ( .NOT. ALLOCATED(lon_full) ) ALLOCATE(lon_full(iim_full,jjm_full))
660	IF ( .NOT. ALLOCATED(lat_full) ) ALLOCATE(lat_full(iim_full,jjm_full))
661	IF ( .NOT. ALLOCATED(mask_full) ) ALLOCATE(mask_full(iim_full,jjm_full))
662	IF ( .NOT. ALLOCATED(mapfac_x_full) ) ALLOCATE(mapfac_x_full(iim_full,jjm_full))
663	IF ( .NOT. ALLOCATED(mapfac_y_full) ) ALLOCATE(mapfac_y_full(iim_full,jjm_full))
664	!
665	DO iv = 1,nvars
666	!
667	iret = NF90_INQUIRE_VARIABLE(fid, iv, name=varname, ndims=varndim, dimids=vardims)
668	!
669	SELECT CASE(varname)
670	!
671	CASE("XLONG_M")
672	iret = NF90_GET_VAR(fid, iv, lon_full)
673	IF (iret /= NF90_NOERR) THEN
674	CALL ipslerr (3,'globgrd_getwrf',"Could not read the longitude from the WRF grid file.", " ", " ")
675	ENDIF
676	CASE("XLAT_M")
677	iret = NF90_GET_VAR(fid, iv, lat_full)
678	IF (iret /= NF90_NOERR) THEN
679	CALL ipslerr (3,'globgrd_getwrf',"Could not read the latitude from the WRF grid file.", " ", " ")
680	ENDIF
681	CASE("LANDMASK")
682	iret = NF90_GET_VAR(fid, iv, mask_full)
683	IF (iret /= NF90_NOERR) THEN
684	CALL ipslerr (3,'globgrd_getwrf',"Could not read the land mask from the WRF grid file.", " ", " ")
685	ENDIF
686	CASE("MAPFAC_MX")
687	iret = NF90_GET_VAR(fid, iv, mapfac_x_full)
688	IF (iret /= NF90_NOERR) THEN
689	CALL ipslerr (3,'globgrd_getwrf',"Could not read the land mask from the WRF grid file.", " ", " ")
690	ENDIF
691	CASE("MAPFAC_MY")
692	iret = NF90_GET_VAR(fid, iv, mapfac_y_full)
693	IF (iret /= NF90_NOERR) THEN
694	CALL ipslerr (3,'globgrd_getwrf',"Could not read the land mask from the WRF grid file.", " ", " ")
695	ENDIF
696	!
697	END SELECT
698	ENDDO
699	!
700	IF ( PRESENT(zoom_lon) .AND. PRESENT(zoom_lat) ) THEN
701	!
702	! if zoomed region
703	!
704	loczoom_lon = zoom_lon
705	loczoom_lat = zoom_lat
706
707
708	! to get new iim, jjm and nbland for the zoomed region
709	! tmp_lon_ind1, tmp_lat_ind1 etc: dimension(2),
710	! with the first one for longitude, second one for latitude
711	tmp_lon_ind1 = MINLOC(ABS(lon_full(:,:)-loczoom_lon(1)))
712	tmp_lat_ind1 = MINLOC(ABS(lat_full(:,:)-loczoom_lat(1)))
713
714	tmp_lon_ind2 = MINLOC(ABS(lon_full(:,:)-loczoom_lon(2)))
715	tmp_lat_ind2 = MINLOC(ABS(lat_full(:,:)-loczoom_lat(2)))
716	!
717	! get indices for the zoomed region
718	! lon1, lat2------- lon2,lat2
719	! \| \|
720	! lon1, lat1------- lon2,lat1
721	!
722	iindex_init = tmp_lon_ind1(1)
723	jindex_init = tmp_lat_ind1(2)
724
725	iindex_end= tmp_lon_ind2(1)
726	jindex_end= tmp_lat_ind2(2)
727
728	! allocate variable values for the zoomed region
729	mask(:,:) = mask_full(iindex_init:iindex_end, jindex_init:jindex_end)
730	lat(:,:) = lat_full(iindex_init:iindex_end, jindex_init:jindex_end)
731	lon(:,:) = lon_full(iindex_init:iindex_end, jindex_init:jindex_end)
732
733	iim_tmp = iindex_end - iindex_init +1
734	jjm_tmp = jindex_end - jindex_init +1
735
736
737	ELSE
738	!
739	! if full grids
740	!
741	iindex_init = 1
742	jindex_init = 1
743
744	! define variable values for full grids
745	lat(:,:) = lat_full(:,:)
746	lon(:,:) = lon_full(:,:)
747	mask(:,:) = mask_full(:,:)
748
749	ENDIF
750	!
751	! Compute corners on the iimxjjm full or zoomed grid
752	DO ip=1,iim
753	DO jp=1,jjm
754	i_orig = iindex_init + ip - 1
755	j_orig = jindex_init + jp - 1
756	! Corners
757	CALL grid_tolola(i_orig+0.5, j_orig+0.5, corners(ip,jp,1,1), corners(ip,jp,1,2))
758	CALL grid_tolola(i_orig+0.5, j_orig-0.5, corners(ip,jp,2,1), corners(ip,jp,2,2))
759	CALL grid_tolola(i_orig-0.5, j_orig-0.5, corners(ip,jp,3,1), corners(ip,jp,3,2))
760	CALL grid_tolola(i_orig-0.5, j_orig+0.5, corners(ip,jp,4,1), corners(ip,jp,4,2))
761	!
762	ENDDO
763	ENDDO
764	!
765	! Compute resolution and area on the gathered, full or zoomed grid
766	!
767	k=0
768	!
769	DO jp=1,jjm
770	DO ip=1,iim
771	!
772	!get the right index of zoomed/full region in the original grids
773	!
774	i_orig = iindex_init + ip - 1
775	j_orig = jindex_init + jp - 1
776	!
777	!
778	! compute area
779	coslat = MAX( COS(lat_full(i_orig,j_orig) * pi/180. ), mincos )
780	dx = ABS(corners(ip,jp,2,1) - corners(ip,jp,1,1)) * pi/180. * R_Earth * coslat
781	dy = ABS(corners(ip,jp,1,2) - corners(ip,jp,3,2)) * pi/180. * R_Earth
782	area(ip,jp) = dx*dy
783	!
784	! compute index and contfrac
785	IF ( mask(ip,jp) > 0.5 ) THEN
786	!
787	! index of the points in the local zoomed grid
788	k=k+1
789	lindex(k) = (jp-1)*iim+ip
790	contfrac(k) = 1.0
791	!
792	ENDIF
793	ENDDO
794	ENDDO
795	!
796	iret = NF90_CLOSE (fid)
797	IF (iret /= NF90_NOERR) THEN
798	CALL ipslerr (3,'globgrd_getwrf',"Error closing the WRF grid file :", " ", " ")
799	ENDIF
800	IF ( ALLOCATED(lon_full) ) DEALLOCATE(lon_full)
801	IF ( ALLOCATED(lat_full) ) DEALLOCATE(lat_full)
802	IF ( ALLOCATED(mask_full) ) DEALLOCATE(mask_full)
803	IF ( ALLOCATED(mapfac_x_full) ) DEALLOCATE(mapfac_x_full)
804	IF ( ALLOCATED(mapfac_y_full) ) DEALLOCATE(mapfac_y_full)
805
806	!
807	END SUBROUTINE globgrd_getwrf
808	!!
809	!! =============================================================================================================================
810	!! SUBROUTINE: globgrd_writegrid
811	!!
812	!>\BRIEF Allows to write the grid to a netDF file for later usage by the glogrid module.
813	!!
814	!! DESCRIPTION: This routine will write a grid description to a netCDF file. mask is on the iimxjjm grid while other
815	!! variables are on the gathered grid.
816	!!
817	!! \n
818	!_ ==============================================================================================================================
819	!!
820	!
821	!
822	SUBROUTINE globgrd_writegrid (gridfilename)
823	!
824	! This routine will write a grid description to a netCDF file. mask is on the iimxjjm grid while other
825	! variables are on the gathered grid.
826	!
827	! ARGUMENTS
828	!
829	CHARACTER(LEN=*), INTENT(in) :: gridfilename
830	!
831	! LOCAL Grid description
832	!
833	INTEGER(i_std) :: iim, jjm, nblindex
834	!
835	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lon, lat
836	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask
837	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: area
838	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:):: corners
839	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: contfrac
840	INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: lindex
841	CHARACTER(LEN=20) :: calendar
842	!
843	! LOCAL netCDF and helping variables
844	!
845	INTEGER(i_std) :: iret, fid, i
846	INTEGER(i_std) :: lonid, latid, landdimid, resid, neighid, maskid, nbcornersid
847	INTEGER(i_std) :: londimid, latdimid, contfracid, resolutionid, neighbourid
848	INTEGER(i_std) :: landindexid, areaid, cornerid
849	!
850	! Get the grid size from the forcing module
851	!
852	CALL forcing_givegridsize (iim, jjm, nblindex)
853	WRITE(,) "Dimension of grid for forcing (iim,jjm,nblindex):", iim,jjm,nblindex
854	!
855	! Allocate fields
856	!
857	ALLOCATE(lon(iim,jjm), lat(iim,jjm))
858	ALLOCATE(mask(iim,jjm))
859	ALLOCATE(area(iim,jjm))
860	ALLOCATE(corners(iim,jjm,4,2))
861	ALLOCATE(lindex(nblindex))
862	ALLOCATE(contfrac(nblindex))
863	!
864	! Get the actual grid from the forcing module
865	!
866	CALL forcing_givegrid(lon, lat, mask, area, corners, lindex, contfrac, calendar)
867	!
868	!
869	iret = NF90_CREATE(gridfilename, NF90_WRITE, fid)
870	IF (iret /= NF90_NOERR) THEN
871	CALL ipslerr (3,'globgrd_writegrid',"Error opening the output file :",gridfilename, " ")
872	ENDIF
873	!
874	! Define dimensions
875	!
876	iret = NF90_DEF_DIM(fid,'lon',iim,londimid)
877	iret = NF90_DEF_DIM(fid,'lat',jjm,latdimid)
878	iret = NF90_DEF_DIM(fid,'nbland',nblindex,landdimid)
879	iret = NF90_DEF_DIM(fid,'nbres',2,resid)
880	iret = NF90_DEF_DIM(fid,'nbcorners',4,nbcornersid)
881	!
882	!
883	! We need to verify here that we have a regulat grid befor deciding if we write lon and lat in 1D or 2D !
884	!
885	!
886	iret = NF90_DEF_VAR(fid,"longitude",NF90_REAL4,londimid,lonid)
887	iret = NF90_PUT_ATT(fid,lonid,'standard_name',"longitude")
888	iret = NF90_PUT_ATT(fid,lonid,'units',"degrees_east")
889	iret = NF90_PUT_ATT(fid,lonid,'valid_min',MINVAL(lon))
890	iret = NF90_PUT_ATT(fid,lonid,'valid_max',MAXVAL(lon))
891	iret = NF90_PUT_ATT(fid,lonid,'long_name',"Longitude")
892	!
893	iret = NF90_DEF_VAR(fid,"latitude",NF90_REAL4,latdimid,latid)
894	iret = NF90_PUT_ATT(fid,latid,'standard_name',"latitude")
895	iret = NF90_PUT_ATT(fid,latid,'units',"degrees_north")
896	iret = NF90_PUT_ATT(fid,latid,'valid_min',MINVAL(lat))
897	iret = NF90_PUT_ATT(fid,latid,'valid_max',MAXVAL(lat))
898	iret = NF90_PUT_ATT(fid,latid,'long_name',"Latitude")
899	!
900	iret = NF90_DEF_VAR(fid,"mask",NF90_REAL4,(/lonid,latid/),maskid)
901	iret = NF90_PUT_ATT(fid,maskid,'standard_name',"mask")
902	iret = NF90_PUT_ATT(fid,maskid,'units',"-")
903	iret = NF90_PUT_ATT(fid,maskid,'valid_min',MINVAL(mask))
904	iret = NF90_PUT_ATT(fid,maskid,'valid_max',MAXVAL(mask))
905	iret = NF90_PUT_ATT(fid,maskid,'long_name',"Land surface mask")
906	!
907	iret = NF90_DEF_VAR(fid,"area",NF90_REAL4,(/lonid,latid/), areaid)
908	iret = NF90_PUT_ATT(fid,areaid,'standard_name',"area")
909	iret = NF90_PUT_ATT(fid,areaid,'units',"m*m")
910	iret = NF90_PUT_ATT(fid,areaid,'valid_min',MINVAL(area))
911	iret = NF90_PUT_ATT(fid,areaid,'valid_max',MAXVAL(area))
912	iret = NF90_PUT_ATT(fid,areaid,'long_name',"Area of grid box")
913	!
914	iret = NF90_DEF_VAR(fid,"corners",NF90_REAL4,(/lonid,latid,nbcornersid,resid/), cornerid)
915	iret = NF90_PUT_ATT(fid,cornerid,'standard_name',"gridcorners")
916	iret = NF90_PUT_ATT(fid,cornerid,'units',"m*m")
917	iret = NF90_PUT_ATT(fid,cornerid,'valid_min',MINVAL(corners))
918	iret = NF90_PUT_ATT(fid,cornerid,'valid_max',MAXVAL(corners))
919	iret = NF90_PUT_ATT(fid,cornerid,'long_name',"corners of grid boxes")
920	!
921	iret = NF90_DEF_VAR(fid,"landindex",NF90_INT, landdimid, landindexid)
922	iret = NF90_PUT_ATT(fid,landindexid,'standard_name',"landindex")
923	iret = NF90_PUT_ATT(fid,landindexid,'units',"-")
924	iret = NF90_PUT_ATT(fid,landindexid,'valid_min',MINVAL(lindex))
925	iret = NF90_PUT_ATT(fid,landindexid,'valid_max',MAXVAL(lindex))
926	iret = NF90_PUT_ATT(fid,landindexid,'long_name',"Land index on global grid (FORTRAN convention)")
927	!
928	iret = NF90_DEF_VAR(fid,"contfrac",NF90_INT,(/landdimid/), contfracid)
929	iret = NF90_PUT_ATT(fid,contfracid,'standard_name',"contfrac")
930	iret = NF90_PUT_ATT(fid,contfracid,'units',"-")
931	iret = NF90_PUT_ATT(fid,contfracid,'valid_min',MINVAL(contfrac))
932	iret = NF90_PUT_ATT(fid,contfracid,'valid_max',MAXVAL(contfrac))
933	iret = NF90_PUT_ATT(fid,contfracid,'long_name',"Fraction of continent in grid box")
934	!
935	! Global attributes
936	!
937	iret = NF90_PUT_ATT(fid, NF90_GLOBAL,'calendar', calendar)
938	!
939	iret = NF90_ENDDEF (fid)
940	IF (iret /= NF90_NOERR) THEN
941	CALL ipslerr (3,'globgrd_writegrid',"Error ending definitions in file :",gridfilename, " ")
942	ENDIF
943	!
944	! Write variables
945	!
946	iret = NF90_PUT_VAR(fid, lonid, lon(:,1))
947	iret = NF90_PUT_VAR(fid, latid, lat(1,:))
948	iret = NF90_PUT_VAR(fid, maskid, mask)
949	iret = NF90_PUT_VAR(fid, areaid, area)
950	iret = NF90_PUT_VAR(fid, cornerid, corners)
951	!
952	iret = NF90_PUT_VAR(fid, landindexid,lindex)
953	iret = NF90_PUT_VAR(fid, contfracid, contfrac)
954	!
955	! Close file
956	!
957	iret = NF90_CLOSE(fid)
958	IF (iret /= NF90_NOERR) THEN
959	CALL ipslerr (3,'globgrd_writegrid',"Error closing the output file :",gridfilename, " ")
960	ENDIF
961	!
962	END SUBROUTINE globgrd_writegrid
963	!!
964	!! =============================================================================================================================
965	!! SUBROUTINE: globgrd_writevar
966	!!
967	!>\BRIEF Writes the grid and a variable to a netCDF file to check if the grid was correctly interpreted by the module.
968	!!
969	!! DESCRIPTION:
970	!!
971	!! \n
972	!_ ==============================================================================================================================
973	!!
974	!
975	SUBROUTINE globgrd_writevar(iim, jjm, lon, lat, nbpt, lalo, var, varname, filename)
976	!
977	! Subroutine used to dump a compressed variable into a full lat/lon grid of a netCDF file
978	!
979	USE netcdf
980	!
981	! ARGUMENTS
982	!
983	INTEGER(i_std), INTENT(in) :: iim, jjm, nbpt
984	REAL(r_std), INTENT(in) :: lon(iim,jjm), lat(iim,jjm)
985	REAL(r_std), INTENT(in) :: lalo(nbpt,2)
986	REAL(r_std), INTENT(in) :: var(nbpt)
987	CHARACTER(LEN=*), INTENT(in) :: varname
988	CHARACTER(LEN=*), INTENT(in) :: filename
989	!
990	! LOCAL
991	!
992	INTEGER(i_std) :: iret, fid, i, ii, jj, nlonid, nlatid, varid
993	INTEGER(i_std) :: ip1, im1, jp1, jm1, di, dj
994	REAL(r_std) :: limlon, limlat
995	INTEGER(i_std), DIMENSION(2) :: lolaid
996	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: varfull, dist
997	INTEGER(i_std), DIMENSION(2) :: closest
998	REAL(r_std), PARAMETER :: epsilon=0.001
999	!
1000	!
1001	WRITE(,) "globgrd_writevar WRITE ", TRIM(varname), " into file ", TRIM(filename)
1002	!
1003	ALLOCATE(varfull(iim,jjm), dist(iim,jjm))
1004	varfull(:,:) = nf90_fill_real
1005	!
1006	! Locate each point on the global grid
1007	!
1008	DO i=1,nbpt
1009	closest(1) = 99999999
1010	closest(2) = 99999999
1011	DO ii=1,iim
1012	DO jj=1,jjm
1013	! Get neighbours
1014	ip1=ii+1
1015	im1=ii-1
1016	jp1=jj+1
1017	jm1=jj-1
1018	di=2
1019	dj=2
1020	! Treat exceptions
1021	IF (ip1 > iim) THEN
1022	ip1=iim
1023	di=1
1024	ENDIF
1025	IF (im1 < 1) THEN
1026	im1=1
1027	di=1
1028	ENDIF
1029	IF ( jp1 > jjm) THEN
1030	jp1=jjm
1031	dj=1
1032	ENDIF
1033	IF ( jm1 < 1) THEN
1034	jm1=1
1035	dj=1
1036	ENDIF
1037	! Calculate limits
1038	limlon=ABS(lon(ip1,jj)-lon(im1,jj))/di-epsilon
1039	limlat=ABS(lat(ii,jp1)-lat(ii,jm1))/dj-epsilon
1040	!
1041	IF ( ABS(lalo(i,1)-lat(ii,jj)) < limlat .AND. ABS(lalo(i,2)-lon(ii,jj)) < limlon ) THEN
1042	closest(1) = ii
1043	closest(2) = jj
1044	ENDIF
1045	ENDDO
1046	ENDDO
1047	IF ( closest(1) > 99999998 .OR. closest(2) > 99999998 ) THEN
1048	WRITE(,) "LALO closest : ", closest
1049	STOP "ERROR in globgrd_writevar"
1050	ELSE
1051	varfull(closest(1),closest(2)) = var(i)
1052	ENDIF
1053	ENDDO
1054	!
1055	! Write the full variable into a NETCDF file
1056	!
1057	iret = NF90_CREATE(filename, NF90_WRITE, fid)
1058	IF (iret /= NF90_NOERR) THEN
1059	CALL ipslerr (3,'globgrd_writevar',"Error opening the output file :",filename, " ")
1060	ENDIF
1061	!
1062	! Define dimensions
1063	!
1064	iret = NF90_DEF_DIM(fid,'Longitude',iim,lolaid(1))
1065	iret = NF90_DEF_DIM(fid,'Latitude',jjm,lolaid(2))
1066	!
1067	iret = NF90_DEF_VAR(fid,"Longitude",NF90_REAL4,lolaid,nlonid)
1068	iret = NF90_PUT_ATT(fid,nlonid,'standard_name',"longitude")
1069	iret = NF90_PUT_ATT(fid,nlonid,'units',"degrees_east")
1070	iret = NF90_PUT_ATT(fid,nlonid,'valid_min',MINVAL(lon))
1071	iret = NF90_PUT_ATT(fid,nlonid,'valid_max',MAXVAL(lon))
1072	iret = NF90_PUT_ATT(fid,nlonid,'long_name',"Longitude")
1073	!
1074	iret = NF90_DEF_VAR(fid,"Latitude",NF90_REAL4,lolaid,nlatid)
1075	iret = NF90_PUT_ATT(fid,nlatid,'standard_name',"latitude")
1076	iret = NF90_PUT_ATT(fid,nlatid,'units',"degrees_north")
1077	iret = NF90_PUT_ATT(fid,nlatid,'valid_min',MINVAL(lat))
1078	iret = NF90_PUT_ATT(fid,nlatid,'valid_max',MAXVAL(lat))
1079	iret = NF90_PUT_ATT(fid,nlatid,'long_name',"Latitude")
1080	!
1081	iret = NF90_DEF_VAR(fid,varname,NF90_REAL4,lolaid,varid)
1082	iret = NF90_PUT_ATT(fid,varid,'_FillValue',NF90_FILL_REAL)
1083	!
1084	iret = NF90_ENDDEF (fid)
1085	IF (iret /= NF90_NOERR) THEN
1086	CALL ipslerr (3,'globgrd_writevar',"Error ending definitions in file :",filename, " ")
1087	ENDIF
1088	!
1089	! Write variables
1090	!
1091	iret = NF90_PUT_VAR(fid,nlonid,lon)
1092	iret = NF90_PUT_VAR(fid,nlatid,lat)
1093	iret = NF90_PUT_VAR(fid,varid,varfull)
1094	!
1095	! Close file
1096	!
1097	iret = NF90_CLOSE(fid)
1098	IF (iret /= NF90_NOERR) THEN
1099	CALL ipslerr (3,'globgrd_writevar',"Error closing the output file :",filename, " ")
1100	ENDIF
1101	!
1102	DEALLOCATE(varfull,dist)
1103	!
1104	WRITE(,) "globgrd_writevar CLOSE file ", TRIM(filename)
1105	!
1106	END SUBROUTINE globgrd_writevar
1107	!
1108	END MODULE globgrd

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

Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/globgrd.f90 @ 8379

Download in other formats: