Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

globgrd.f90 @ 7541

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

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

Context Navigation

source: branches/publications/ORCHIDEE_2.2_r7266/ORCHIDEE/src_driver/globgrd.f90 @ 7541

Download in other formats: