Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-57/src_driver/globgrd.f90 @ 7746

Last change on this file since 7746 was 3447, checked in by josefine.ghattas, 9 years ago
Integrated branch ORCHIDEE-DRIVER in the trunk. See #244
File size: 31.6 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) :: iim_full, jjm_full, nbland_full
83	CHARACTER(LEN=20) :: axname, varname
84	CHARACTER(LEN=120) :: tmpfile
85	REAL(r_std), DIMENSION(2) :: loczoom_lon, loczoom_lat
86	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask
87	!
88	! Set default values against which we can test
89	!
90	iim = -1
91	jjm = -1
92	!
93	! Verify the grid file name
94	!
95	IF ( INDEX(filename,"NONE") >= 1 ) THEN
96	is_forcing_file=.TRUE.
97	IF ( PRESENT(forcingfile) ) THEN
98	tmpfile=forcingfile(1)
99	ELSE
100	CALL ipslerr (3,'globgrd_getdomsz',"Error No grid file provided :",tmpfile, " ")
101	ENDIF
102	ELSE
103	is_forcing_file=.FALSE.
104	tmpfile=filename
105	ENDIF
106	!
107	! Verify that the zomm is provided. Else choose the entire globe
108	!
109	IF ( PRESENT(zoom_lon) .AND. PRESENT(zoom_lat) ) THEN
110	loczoom_lon = zoom_lon
111	loczoom_lat = zoom_lat
112	ELSE
113	loczoom_lon(1) = -180.0
114	loczoom_lon(2) = 180.0
115	loczoom_lat(1) = -90.0
116	loczoom_lat(2) = 90.0
117	ENDIF
118	!
119	! Open the correct file
120	!
121	iret = NF90_OPEN (tmpfile, NF90_NOWRITE, fid)
122	IF (iret /= NF90_NOERR) THEN
123	CALL ipslerr (3,'globgrd_getdomsz',"Error opening the grid file :",tmpfile, " ")
124	ENDIF
125	!
126	!
127	iret = NF90_INQUIRE (fid, nDimensions=ndims, nVariables=nvars, &
128	nAttributes=nb_atts, unlimitedDimId=id_unlim)
129	IF (iret /= NF90_NOERR) THEN
130	CALL ipslerr (3,'globgrd_getdomsz',"Error in NF90_INQUIRE :",tmpfile, " ")
131	ENDIF
132	!
133	DO iv=1,ndims
134	!
135	iret = NF90_INQUIRE_DIMENSION (fid, iv, name=axname, len=lll)
136	IF (iret /= NF90_NOERR) THEN
137	CALL ipslerr (3,'globgrd_getdomsz',"Could not get size of dimension :"," "," ")
138	ENDIF
139	!
140	! This can be refined by testing the actual grid found in the file.
141	!
142	SELECT CASE(axname)
143	!
144	!! Coordinate variables used by WRF.
145	CASE("west_east")
146	iim = lll
147	model_guess = "WRF"
148	CASE("south_north")
149	jjm = lll
150	model_guess = "WRF"
151	!
152	!! Variables used in WFDEI
153	CASE("lon")
154	iim = lll
155	model_guess = "regular"
156	CASE("lat")
157	jjm = lll
158	model_guess = "regular"
159	CASE("nbland")
160	nbland = lll
161	!
162	!! Variables used by CRU-NCEP
163	CASE("nav_lon")
164	iim = lll
165	model_guess = "regular"
166	CASE("nav_lat")
167	jjm = lll
168	model_guess = "regular"
169	CASE("land")
170	nbland = lll
171	END SELECT
172	ENDDO
173	!
174	! If we have a WRF file we need to count the number of land points.
175	!
176	IF ( model_guess == "WRF" ) THEN
177
178	ALLOCATE(mask(iim,jjm))
179
180	varname = "LANDMASK"
181	iret = NF90_INQ_VARID (fid, varname, iv)
182	IF (iret /= NF90_NOERR) THEN
183	CALL ipslerr (3,'globgrd_getdomsz',"Could not find variable ", varname," ")
184	ELSE
185	iret = NF90_GET_VAR (fid,iv,mask)
186	ENDIF
187
188	nbland = COUNT(mask > 0.5)
189
190	ENDIF
191	!
192	!
193	! If we are in the case of a forcing file a few functions from forcing_tools need to be called
194	! so that the file is analysed with the tools of the forcing module.
195	!
196	IF ( is_forcing_file ) THEN
197	!
198	! Because we are re-using routines from the forcing module, we have to
199	! close the file. It will be opened again by the forcing module.
200	!
201	iret = NF90_CLOSE(fid)
202	IF (iret /= NF90_NOERR) THEN
203	CALL ipslerr (3,'globgrd_getdomzz',"Error closing the output file :",filename, " ")
204	ENDIF
205	!
206	! This has to be a regular grid. A more clever clasification of files will be needed.
207	model_guess = "regular"
208	!
209	CALL forcing_getglogrid (1, forcingfile, iim_full, jjm_full, nbland_full, .TRUE.)
210	WRITE(,) forcingfile, "Forcing file with dimensions : ", iim_full, jjm_full, nbland_full
211	!
212	CALL forcing_zoomgrid (loczoom_lon, loczoom_lat, forcingfile(1), .TRUE.)
213	!
214	CALL forcing_givegridsize (iim, jjm, nbland)
215	!
216	ENDIF
217	!
218	! Do a final test to see if we got the information needed.
219	!
220	IF ( iim < 0 .OR. jjm < 0 ) THEN
221	CALL ipslerr (3,'globgrd_getdomsz',"Could not get the horizontal size of the domaine out of the file",&
222	& filename,"Are you sure that the case for this type of file is foreseen ? ")
223	ENDIF
224	!
225	!
226	END SUBROUTINE globgrd_getdomsz
227	!!
228	!! =============================================================================================================================
229	!! SUBROUTINE: globgrd_getgrid
230	!!
231	!>\BRIEF This routine extracts the coordinates and land/sea mask from the domain files.
232	!!
233	!! DESCRIPTION: The domain size is provided together with the netCDF file ID so that the main information can
234	!! be extracted from the file. We will read the longitude, latitude, land/sea mask and calendar.
235	!! This allows to set-up ORCHIDEE. We also provide the corners of the grid-boxes as this is needed
236	!! for setting-up OASIS but is computed more correctly in grid.f90 for ORCHIDEE.
237	!! This routine is only an interface to globgrd_getwrf, globgrd_getregular and forcing_givegrid.
238	!! forcing_givegrid is an interface to the forcing_tools.f90 module so that we are certain to have
239	!! the same grid information between both modules.
240	!!
241	!! \n
242	!_ ==============================================================================================================================
243	!!
244	!---------------------------------------------------------------------
245	!-
246	!-
247	!---------------------------------------------------------------------
248	SUBROUTINE globgrd_getgrid(fid, iim, jjm, nbland, model_guess, lon, lat, mask, area, corners, &
249	& lindex, contfrac, calendar)
250	!
251	!
252	! This subroutine only switched between routines to extract and compte the grid data needed for
253	! ORCHIDEE.
254	!
255	!
256	! INPUT
257	!
258	INTEGER(i_std), INTENT(in) :: fid
259	INTEGER(i_std), INTENT(in) :: iim, jjm, nbland
260	CHARACTER(LEN=*), INTENT(in) :: model_guess
261	!
262	! OUTPUT
263	!
264	REAL(r_std),DIMENSION(iim,jjm), INTENT(out) :: lon, lat, mask, area
265	REAL(r_std),DIMENSION(iim,jjm,4,2), INTENT(out) :: corners
266	INTEGER(i_std), DIMENSION(nbland), INTENT(out) :: lindex
267	REAL(r_std),DIMENSION(nbland), INTENT(out) :: contfrac
268	CHARACTER(LEN=20), INTENT(out) :: calendar
269	!
270	SELECT CASE(model_guess)
271
272	CASE("WRF")
273	CALL globgrd_getwrf(fid, iim, jjm, nbland, lon, lat, mask, area, corners, &
274	& lindex, contfrac, calendar)
275	CASE("regular")
276	IF ( .NOT. is_forcing_file ) THEN
277	CALL globgrd_getregular(fid, iim, jjm, nbland, lon, lat, mask, area, corners, &
278	& lindex, contfrac, calendar)
279	ELSE
280	CALL forcing_givegrid(lon, lat, mask, area, corners, lindex, contfrac, calendar)
281	CALL forcing_close()
282	ENDIF
283	CASE DEFAULT
284	CALL ipslerr (3,'globgrd_getgrid',"The model/grid type we guessed is not recognized here. model_guess =",&
285	& model_guess,"Have you used the right file and are you sure that this case is foreseen ? ")
286	END SELECT
287	!
288	!
289	END SUBROUTINE globgrd_getgrid
290	!!
291	!! =============================================================================================================================
292	!! SUBROUTINE: globgrd_regular
293	!!
294	!>\BRIEF The routine to obtain regular grids from the file.
295	!!
296	!! DESCRIPTION: Read the regular grid and its information from the opened file (fid).
297	!!
298	!! \n
299	!_ ==============================================================================================================================
300	!!
301	SUBROUTINE globgrd_getregular(fid, iim, jjm, nbland, lon, lat, mask, area, corners, &
302	& lindex, contfrac, calendar)
303	!
304	USE defprec
305	USE netcdf
306	!
307	! INPUT
308	!
309	INTEGER(i_std), INTENT(in) :: fid
310	INTEGER(i_std), INTENT(in) :: iim, jjm, nbland
311	!
312	! OUTPUT
313	!
314	REAL(r_std),DIMENSION(iim,jjm), INTENT(out) :: lon, lat, mask, area
315	REAL(r_std),DIMENSION(iim,jjm,4,2), INTENT(out) :: corners
316	INTEGER(i_std), DIMENSION(nbland), INTENT(out) :: lindex
317	REAL(r_std),DIMENSION(nbland), INTENT(out) :: contfrac
318	CHARACTER(LEN=20), INTENT(out) :: calendar
319	!
320	! LOCAL
321	!
322	INTEGER(i_std) :: iret, iv, nvars, varndim
323	INTEGER(i_std) :: i, j
324	CHARACTER(LEN=20) :: varname
325	INTEGER(i_std), DIMENSION(4) :: vardims
326	REAL(r_std) :: dx
327	!
328	! Set some default values agains which we can check
329	!
330	lon(:,:) = val_exp
331	lat(:,:) = val_exp
332	mask(:,:) = val_exp
333	area(:,:) = val_exp
334	corners(:,:,:,:) = val_exp
335	!
336	lindex(:) = INT(val_exp)
337	contfrac(:) = val_exp
338	!
339	! Get the global attributes from grid file
340	!
341	iret = NF90_GET_ATT(fid, NF90_GLOBAL, 'calendar', calendar)
342	IF (iret /= NF90_NOERR) THEN
343	CALL ipslerr (3,'globgrd_getregular',"Could not read the calendar in grid file.", " ", " ")
344	ENDIF
345	!
346	iret = NF90_INQUIRE (fid, nVariables=nvars)
347	!
348	DO iv = 1,nvars
349	!
350	iret = NF90_INQUIRE_VARIABLE(fid, iv, name=varname, ndims=varndim, dimids=vardims)
351	!
352	!
353	SELECT CASE(varname)
354	!
355	CASE("longitude")
356	IF (varndim == 1 ) THEN
357	DO j=1,jjm
358	iret = NF90_GET_VAR(fid, iv, lon(:,j))
359	ENDDO
360	ELSE IF (varndim == 2 ) THEN
361	iret = NF90_GET_VAR(fid, iv, lon)
362	ELSE
363	CALL ipslerr (3,'globgrd_getregular',"Longitude cannot have more than 2 dimensions","","")
364	ENDIF
365
366	CASE ("latitude")
367	IF (varndim == 1 ) THEN
368	DO i=1,iim
369	iret = NF90_GET_VAR(fid, iv, lat(i,:))
370	ENDDO
371	ELSE IF (varndim == 2 ) THEN
372	iret = NF90_GET_VAR(fid, iv, lon)
373	ELSE
374	CALL ipslerr (3,'globgrd_getregular',"Latitude cannot have more than 2 dimensions","","")
375	ENDIF
376
377	CASE ("mask")
378	IF (varndim /= 2 ) THEN
379	CALL ipslerr (3,'globgrd_getregular',"mask needs to have 2 dimensions","","")
380	ELSE
381	iret = NF90_GET_VAR (fid, iv, mask)
382	ENDIF
383
384	CASE ("areas")
385	IF (varndim /= 2 ) THEN
386	CALL ipslerr (3,'globgrd_getregular',"Areas needs to have 2 dimensions","","")
387	ELSE
388	iret = NF90_GET_VAR (fid, iv, area)
389	ENDIF
390
391	CASE ("corners")
392	IF (varndim /= 4 ) THEN
393	CALL ipslerr (3,'globgrd_getregular',"corners needs to have 4 dimensions","","")
394	ELSE
395	iret = NF90_GET_VAR (fid, iv, corners)
396	ENDIF
397
398	CASE ("landindex")
399	IF (varndim /= 1 ) THEN
400	CALL ipslerr (3,'globgrd_getregular',"landindex is the list of continental points to be gathered", &
401	& "Thus it can only have 1 dimensions","")
402	ELSE
403	iret = NF90_GET_VAR (fid, iv, lindex)
404	ENDIF
405
406	CASE ("contfrac")
407	IF (varndim /= 1 ) THEN
408	CALL ipslerr (3,'globgrd_getregular',"Contfrac needs to be a gathered variable", &
409	& "thus it needs only 1 dimensions","")
410	ELSE
411	iret = NF90_GET_VAR (fid, iv, contfrac)
412	ENDIF
413
414	END SELECT
415	!
416	ENDDO
417	!
418	!
419	iret = NF90_CLOSE(fid)
420	!
421	! Verify that we have al the variables needed to describe the ORCHIDEE grid
422	!
423	IF ( ANY( lon(:,:) == val_exp ) ) THEN
424	CALL ipslerr (3,'globgrd_getregular',"The longitude of the ORCHIDEE grid could not be extracted from the",&
425	& "grid definition file","")
426	ENDIF
427	!
428	IF ( ANY( lat(:,:) == val_exp ) ) THEN
429	CALL ipslerr (3,'globgrd_getregular',"The latitude of the ORCHIDEE grid could not be extracted from the",&
430	& "grid definition file","")
431	ENDIF
432	!
433	IF ( ANY( lindex(:) == INT(val_exp) ) ) THEN
434	CALL ipslerr (3,'globgrd_getregular',"The lindex of the ORCHIDEE grid could not be extracted from the",&
435	& "grid definition file","")
436	ENDIF
437	!
438	IF ( ALL( mask(:,:) == val_exp ) ) THEN
439	CALL ipslerr (3,'globgrd_getregular',"The land mask of the ORCHIDEE grid could not be extracted from the",&
440	& "grid definition file","")
441	ELSE IF (MAXVAL(mask) > 1 ) THEN
442	CALL ipslerr (2,'globgrd_getregular',"We have a special case for the mask which needs to be treated.",&
443	& "The field contains the indices of the land points on a compressed grid.","So we replace them with 1 or 0.")
444	DO i=1,iim
445	DO j=1,jjm
446	IF ( mask(i,j) > iim*jjm ) THEN
447	mask(i,j) = 0
448	ELSE
449	mask(i,j) = 1
450	ENDIF
451	ENDDO
452	ENDDO
453	ENDIF
454	!
455	IF ( ANY( contfrac(:) == val_exp ) ) THEN
456	CALL ipslerr (2,'globgrd_getregular',"The continental fraction of the ORCHIDEE grid could not be extracted from the",&
457	& "grid definition file","Thus on all land points it is set to 1.")
458	contfrac(:) = 1.
459	ENDIF
460	!
461	IF ( ANY( corners(:,:,:,:) == val_exp ) ) THEN
462	CALL ipslerr (3,'globgrd_getregular',"The corners for the ORCHIDEE grid could not be extracted from the",&
463	& "grid definition file","As we have to assume a very general grid we cannot do anything !")
464	ENDIF
465	!
466	!
467	END SUBROUTINE globgrd_getregular
468	!!
469	!! =============================================================================================================================
470	!! SUBROUTINE: globgrd_getwrf
471	!!
472	!>\BRIEF Routine to read the WRF grid description file.
473	!!
474	!! DESCRIPTION: Read the WRF grid and its information from the opened file (fid) and convert
475	!! it to the variables needed by ORCHIDEE.
476	!!
477	!! \n
478	!_ ==============================================================================================================================
479	!!
480	SUBROUTINE globgrd_getwrf(fid, iim, jjm, nbland, lon, lat, mask, area, corners, &
481	& lindex, contfrac, calendar)
482	!
483	USE defprec
484	USE netcdf
485	!
486	! INPUT
487	!
488	INTEGER(i_std), INTENT(in) :: fid
489	INTEGER(i_std), INTENT(in) :: iim, jjm, nbland
490	!
491	! OUTPUT
492	!
493	REAL(r_std),DIMENSION(iim,jjm), INTENT(out) :: lon, lat, mask, area
494	REAL(r_std),DIMENSION(iim,jjm,4,2), INTENT(out) :: corners
495	INTEGER(i_std), DIMENSION(nbland), INTENT(out) :: lindex
496	REAL(r_std),DIMENSION(nbland), INTENT(out) :: contfrac
497	CHARACTER(LEN=20), INTENT(out) :: calendar
498	!
499	! LOCAL
500	!
501	INTEGER(i_std) :: i, ip, jp, k, iret, iv, nvars, varndim
502	INTEGER(i_std) :: imm1, imp1
503	CHARACTER(LEN=20) :: varname
504	REAL(r_std),DIMENSION(iim,jjm) :: mapfac_x, mapfac_y
505	INTEGER(i_std), DIMENSION(4) :: vardims
506	INTEGER(i_std), DIMENSION(8) :: rose
507	!
508	!
509	! Set some default values agains which we can check
510	!
511	lon(:,:) = val_exp
512	lat(:,:) = val_exp
513	mask(:,:) = val_exp
514	area(:,:) = val_exp
515	corners(:,:,:,:) = val_exp
516	!
517	lindex(:) = INT(val_exp)
518	contfrac(:) = val_exp
519	!
520	calendar = 'gregorian'
521	!
522	!
523	! Init projection in grid.f90 so that it can be used later for projections.
524	!
525	CALL grid_initproj(fid, iim, jjm)
526	!
527	iret = NF90_INQUIRE (fid, nVariables=nvars)
528	!
529	IF (iret /= NF90_NOERR) THEN
530	CALL ipslerr (3,'globgrd_getwrf',"Error inquiering variables from WRF grid file."," ", " ")
531	ENDIF
532	!
533	DO iv = 1,nvars
534	!
535	iret = NF90_INQUIRE_VARIABLE(fid, iv, name=varname, ndims=varndim, dimids=vardims)
536	!
537	SELECT CASE(varname)
538	!
539	CASE("XLONG_M")
540	iret = NF90_GET_VAR(fid, iv, lon)
541	IF (iret /= NF90_NOERR) THEN
542	CALL ipslerr (3,'globgrd_getwrf',"Could not read the longitude from the WRF grid file.", " ", " ")
543	ENDIF
544	CASE("XLAT_M")
545	iret = NF90_GET_VAR(fid, iv, lat)
546	IF (iret /= NF90_NOERR) THEN
547	CALL ipslerr (3,'globgrd_getwrf',"Could not read the latitude from the WRF grid file.", " ", " ")
548	ENDIF
549	CASE("LANDMASK")
550	iret = NF90_GET_VAR(fid, iv, mask)
551	IF (iret /= NF90_NOERR) THEN
552	CALL ipslerr (3,'globgrd_getwrf',"Could not read the land mask from the WRF grid file.", " ", " ")
553	ENDIF
554	CASE("MAPFAC_MX")
555	iret = NF90_GET_VAR(fid, iv, mapfac_x)
556	IF (iret /= NF90_NOERR) THEN
557	CALL ipslerr (3,'globgrd_getwrf',"Could not read the land mask from the WRF grid file.", " ", " ")
558	ENDIF
559	CASE("MAPFAC_MY")
560	iret = NF90_GET_VAR(fid, iv, mapfac_y)
561	IF (iret /= NF90_NOERR) THEN
562	CALL ipslerr (3,'globgrd_getwrf',"Could not read the land mask from the WRF grid file.", " ", " ")
563	ENDIF
564	!
565	END SELECT
566	ENDDO
567	!
568	! Compute corners and area on the iimxjjm grid
569	!
570	DO ip=1,iim
571	DO jp=1,jjm
572	! Corners
573	CALL grid_tolola(ip+0.5, jp+0.5, corners(ip,jp,1,1), corners(ip,jp,1,2))
574	CALL grid_tolola(ip+0.5, jp-0.5, corners(ip,jp,2,1), corners(ip,jp,2,2))
575	CALL grid_tolola(ip-0.5, jp-0.5, corners(ip,jp,3,1), corners(ip,jp,3,2))
576	CALL grid_tolola(ip-0.5, jp-0.5, corners(ip,jp,4,1), corners(ip,jp,4,2))
577	!
578	ENDDO
579	ENDDO
580	!
581	! Compute resolution and area on the gathered grid
582	!
583	k=0
584	!
585	DO jp=1,jjm
586	DO ip=1,iim
587	IF ( mask(ip,jp) > 0.5 ) THEN
588	!
589	k=k+1
590	lindex(k) = (jp-1)*iim+ip
591	contfrac(k) = 1.0
592	!
593	ENDIF
594	ENDDO
595	ENDDO
596	!
597	iret = NF90_CLOSE (fid)
598	IF (iret /= NF90_NOERR) THEN
599	CALL ipslerr (3,'globgrd_getwrf',"Error closing the WRF grid file :", " ", " ")
600	ENDIF
601	!
602	END SUBROUTINE globgrd_getwrf
603	!!
604	!! =============================================================================================================================
605	!! SUBROUTINE: globgrd_writegrid
606	!!
607	!>\BRIEF Allows to write the grid to a netDF file for later usage by the glogrid module.
608	!!
609	!! DESCRIPTION: This routine will write a grid description to a netCDF file. mask is on the iimxjjm grid while other
610	!! variables are on the gathered grid.
611	!!
612	!! \n
613	!_ ==============================================================================================================================
614	!!
615	!
616	!
617	SUBROUTINE globgrd_writegrid (gridfilename)
618	!
619	! This routine will write a grid description to a netCDF file. mask is on the iimxjjm grid while other
620	! variables are on the gathered grid.
621	!
622	! ARGUMENTS
623	!
624	CHARACTER(LEN=*), INTENT(in) :: gridfilename
625	!
626	! LOCAL Grid description
627	!
628	INTEGER(i_std) :: iim, jjm, nblindex
629	!
630	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lon, lat
631	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask
632	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: area
633	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:):: corners
634	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: contfrac
635	INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: lindex
636	CHARACTER(LEN=20) :: calendar
637	!
638	! LOCAL netCDF and helping variables
639	!
640	INTEGER(i_std) :: iret, fid, i
641	INTEGER(i_std) :: lonid, latid, landdimid, resid, neighid, maskid, nbcornersid
642	INTEGER(i_std) :: londimid, latdimid, contfracid, resolutionid, neighbourid
643	INTEGER(i_std) :: landindexid, areaid, cornerid
644	!
645	! Get the grid size from the forcing module
646	!
647	CALL forcing_givegridsize (iim, jjm, nblindex)
648	WRITE(,) "Dimension of grid for forcing (iim,jjm,nblindex):", iim,jjm,nblindex
649	!
650	! Allocate fields
651	!
652	ALLOCATE(lon(iim,jjm), lat(iim,jjm))
653	ALLOCATE(mask(iim,jjm))
654	ALLOCATE(area(iim,jjm))
655	ALLOCATE(corners(iim,jjm,4,2))
656	ALLOCATE(lindex(nblindex))
657	ALLOCATE(contfrac(nblindex))
658	!
659	! Get the actual grid from the forcing module
660	!
661	CALL forcing_givegrid(lon, lat, mask, area, corners, lindex, contfrac, calendar)
662	!
663	!
664	iret = NF90_CREATE(gridfilename, NF90_WRITE, fid)
665	IF (iret /= NF90_NOERR) THEN
666	CALL ipslerr (3,'globgrd_writegrid',"Error opening the output file :",gridfilename, " ")
667	ENDIF
668	!
669	! Define dimensions
670	!
671	iret = NF90_DEF_DIM(fid,'lon',iim,londimid)
672	iret = NF90_DEF_DIM(fid,'lat',jjm,latdimid)
673	iret = NF90_DEF_DIM(fid,'nbland',nblindex,landdimid)
674	iret = NF90_DEF_DIM(fid,'nbres',2,resid)
675	iret = NF90_DEF_DIM(fid,'nbcorners',4,nbcornersid)
676	!
677	!
678	! We need to verify here that we have a regulat grid befor deciding if we write lon and lat in 1D or 2D !
679	!
680	!
681	iret = NF90_DEF_VAR(fid,"longitude",NF90_REAL4,londimid,lonid)
682	iret = NF90_PUT_ATT(fid,lonid,'standard_name',"longitude")
683	iret = NF90_PUT_ATT(fid,lonid,'units',"degrees_east")
684	iret = NF90_PUT_ATT(fid,lonid,'valid_min',MINVAL(lon))
685	iret = NF90_PUT_ATT(fid,lonid,'valid_max',MAXVAL(lon))
686	iret = NF90_PUT_ATT(fid,lonid,'long_name',"Longitude")
687	!
688	iret = NF90_DEF_VAR(fid,"latitude",NF90_REAL4,latdimid,latid)
689	iret = NF90_PUT_ATT(fid,latid,'standard_name',"latitude")
690	iret = NF90_PUT_ATT(fid,latid,'units',"degrees_north")
691	iret = NF90_PUT_ATT(fid,latid,'valid_min',MINVAL(lat))
692	iret = NF90_PUT_ATT(fid,latid,'valid_max',MAXVAL(lat))
693	iret = NF90_PUT_ATT(fid,latid,'long_name',"Latitude")
694	!
695	iret = NF90_DEF_VAR(fid,"mask",NF90_REAL4,(/lonid,latid/),maskid)
696	iret = NF90_PUT_ATT(fid,maskid,'standard_name',"mask")
697	iret = NF90_PUT_ATT(fid,maskid,'units',"-")
698	iret = NF90_PUT_ATT(fid,maskid,'valid_min',MINVAL(mask))
699	iret = NF90_PUT_ATT(fid,maskid,'valid_max',MAXVAL(mask))
700	iret = NF90_PUT_ATT(fid,maskid,'long_name',"Land surface mask")
701	!
702	iret = NF90_DEF_VAR(fid,"area",NF90_REAL4,(/lonid,latid/), areaid)
703	iret = NF90_PUT_ATT(fid,areaid,'standard_name',"area")
704	iret = NF90_PUT_ATT(fid,areaid,'units',"m*m")
705	iret = NF90_PUT_ATT(fid,areaid,'valid_min',MINVAL(area))
706	iret = NF90_PUT_ATT(fid,areaid,'valid_max',MAXVAL(area))
707	iret = NF90_PUT_ATT(fid,areaid,'long_name',"Area of grid box")
708	!
709	iret = NF90_DEF_VAR(fid,"corners",NF90_REAL4,(/lonid,latid,nbcornersid,resid/), cornerid)
710	iret = NF90_PUT_ATT(fid,cornerid,'standard_name',"gridcorners")
711	iret = NF90_PUT_ATT(fid,cornerid,'units',"m*m")
712	iret = NF90_PUT_ATT(fid,cornerid,'valid_min',MINVAL(corners))
713	iret = NF90_PUT_ATT(fid,cornerid,'valid_max',MAXVAL(corners))
714	iret = NF90_PUT_ATT(fid,cornerid,'long_name',"corners of grid boxes")
715	!
716	iret = NF90_DEF_VAR(fid,"landindex",NF90_INT, landdimid, landindexid)
717	iret = NF90_PUT_ATT(fid,landindexid,'standard_name',"landindex")
718	iret = NF90_PUT_ATT(fid,landindexid,'units',"-")
719	iret = NF90_PUT_ATT(fid,landindexid,'valid_min',MINVAL(lindex))
720	iret = NF90_PUT_ATT(fid,landindexid,'valid_max',MAXVAL(lindex))
721	iret = NF90_PUT_ATT(fid,landindexid,'long_name',"Land index on global grid (FORTRAN convention)")
722	!
723	iret = NF90_DEF_VAR(fid,"contfrac",NF90_INT,(/landdimid/), contfracid)
724	iret = NF90_PUT_ATT(fid,contfracid,'standard_name',"contfrac")
725	iret = NF90_PUT_ATT(fid,contfracid,'units',"-")
726	iret = NF90_PUT_ATT(fid,contfracid,'valid_min',MINVAL(contfrac))
727	iret = NF90_PUT_ATT(fid,contfracid,'valid_max',MAXVAL(contfrac))
728	iret = NF90_PUT_ATT(fid,contfracid,'long_name',"Fraction of continent in grid box")
729	!
730	! Global attributes
731	!
732	iret = NF90_PUT_ATT(fid, NF90_GLOBAL,'calendar', calendar)
733	!
734	iret = NF90_ENDDEF (fid)
735	IF (iret /= NF90_NOERR) THEN
736	CALL ipslerr (3,'globgrd_writegrid',"Error ending definitions in file :",gridfilename, " ")
737	ENDIF
738	!
739	! Write variables
740	!
741	iret = NF90_PUT_VAR(fid, lonid, lon(:,1))
742	iret = NF90_PUT_VAR(fid, latid, lat(1,:))
743	iret = NF90_PUT_VAR(fid, maskid, mask)
744	iret = NF90_PUT_VAR(fid, areaid, area)
745	iret = NF90_PUT_VAR(fid, cornerid, corners)
746	!
747	iret = NF90_PUT_VAR(fid, landindexid,lindex)
748	iret = NF90_PUT_VAR(fid, contfracid, contfrac)
749	!
750	! Close file
751	!
752	iret = NF90_CLOSE(fid)
753	IF (iret /= NF90_NOERR) THEN
754	CALL ipslerr (3,'globgrd_writegrid',"Error closing the output file :",gridfilename, " ")
755	ENDIF
756	!
757	END SUBROUTINE globgrd_writegrid
758	!!
759	!! =============================================================================================================================
760	!! SUBROUTINE: globgrd_writevar
761	!!
762	!>\BRIEF Writes the grid and a variable to a netCDF file to check if the grid was correctly interpreted by the module.
763	!!
764	!! DESCRIPTION:
765	!!
766	!! \n
767	!_ ==============================================================================================================================
768	!!
769	!
770	SUBROUTINE globgrd_writevar(iim, jjm, lon, lat, nbpt, lalo, var, varname, filename)
771	!
772	! Subroutine used to dump a compressed variable into a full lat/lon grid of a netCDF file
773	!
774	USE netcdf
775	!
776	! ARGUMENTS
777	!
778	INTEGER(i_std), INTENT(in) :: iim, jjm, nbpt
779	REAL(r_std), INTENT(in) :: lon(iim,jjm), lat(iim,jjm)
780	REAL(r_std), INTENT(in) :: lalo(nbpt,2)
781	REAL(r_std), INTENT(in) :: var(nbpt)
782	CHARACTER(LEN=*), INTENT(in) :: varname
783	CHARACTER(LEN=*), INTENT(in) :: filename
784	!
785	! LOCAL
786	!
787	INTEGER(i_std) :: iret, fid, i, ii, jj, nlonid, nlatid, varid
788	INTEGER(i_std), DIMENSION(2) :: lolaid
789	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: varfull, dist
790	INTEGER(i_std), DIMENSION(2) :: closest
791	!
792	!
793	WRITE(,) "globgrd_writevar WRITE ", TRIM(varname), " into file ", TRIM(filename)
794	!
795	ALLOCATE(varfull(iim,jjm), dist(iim,jjm))
796	varfull(:,:) = nf90_fill_real
797	!
798	! Locate each point on the global grid
799	!
800	DO i=1,nbpt
801	closest(1) = 99999999
802	closest(2) = 99999999
803	DO ii=1,iim
804	DO jj=1,jjm
805	IF ( ABS(lalo(i,1)-lat(ii,jj)) < 0.25 .AND. ABS(lalo(i,2)-lon(ii,jj)) < 0.25 ) THEN
806	closest(1) = ii
807	closest(2) = jj
808	ENDIF
809	ENDDO
810	ENDDO
811	IF ( closest(1) > 99999998 .OR. closest(2) > 99999998 ) THEN
812	WRITE(,) "LALO closest : ", closest
813	STOP "ERROR in globgrd_writevar"
814	ELSE
815	varfull(closest(1),closest(2)) = var(i)
816	ENDIF
817	ENDDO
818	!
819	! Write the full variable into a NETCDF file
820	!
821	iret = NF90_CREATE(filename, NF90_WRITE, fid)
822	IF (iret /= NF90_NOERR) THEN
823	CALL ipslerr (3,'globgrd_writevar',"Error opening the output file :",filename, " ")
824	ENDIF
825	!
826	! Define dimensions
827	!
828	iret = NF90_DEF_DIM(fid,'Longitude',iim,lolaid(1))
829	iret = NF90_DEF_DIM(fid,'Latitude',jjm,lolaid(2))
830	!
831	iret = NF90_DEF_VAR(fid,"Longitude",NF90_REAL4,lolaid,nlonid)
832	iret = NF90_PUT_ATT(fid,nlonid,'standard_name',"longitude")
833	iret = NF90_PUT_ATT(fid,nlonid,'units',"degrees_east")
834	iret = NF90_PUT_ATT(fid,nlonid,'valid_min',MINVAL(lon))
835	iret = NF90_PUT_ATT(fid,nlonid,'valid_max',MAXVAL(lon))
836	iret = NF90_PUT_ATT(fid,nlonid,'long_name',"Longitude")
837	!
838	iret = NF90_DEF_VAR(fid,"Latitude",NF90_REAL4,lolaid,nlatid)
839	iret = NF90_PUT_ATT(fid,nlatid,'standard_name',"latitude")
840	iret = NF90_PUT_ATT(fid,nlatid,'units',"degrees_north")
841	iret = NF90_PUT_ATT(fid,nlatid,'valid_min',MINVAL(lat))
842	iret = NF90_PUT_ATT(fid,nlatid,'valid_max',MAXVAL(lat))
843	iret = NF90_PUT_ATT(fid,nlatid,'long_name',"Latitude")
844	!
845	iret = NF90_DEF_VAR(fid,varname,NF90_REAL4,lolaid,varid)
846	iret = NF90_PUT_ATT(fid,varid,'_FillValue',NF90_FILL_REAL)
847	!
848	iret = NF90_ENDDEF (fid)
849	IF (iret /= NF90_NOERR) THEN
850	CALL ipslerr (3,'globgrd_writevar',"Error ending definitions in file :",filename, " ")
851	ENDIF
852	!
853	! Write variables
854	!
855	iret = NF90_PUT_VAR(fid,nlonid,lon)
856	iret = NF90_PUT_VAR(fid,nlatid,lat)
857	iret = NF90_PUT_VAR(fid,varid,varfull)
858	!
859	! Close file
860	!
861	iret = NF90_CLOSE(fid)
862	IF (iret /= NF90_NOERR) THEN
863	CALL ipslerr (3,'globgrd_writevar',"Error closing the output file :",filename, " ")
864	ENDIF
865	!
866	DEALLOCATE(varfull,dist)
867	!
868	WRITE(,) "globgrd_writevar CLOSE file ", TRIM(filename)
869	!
870	END SUBROUTINE globgrd_writevar
871	!
872	END MODULE globgrd

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

Download in other formats: