Context Navigation

source: branches/publications/ORCHIDEE_IPSLCM5A2.1.r5307/ORCHIDEE/src_global/grid.f90 @ 7683

Last change on this file since 7683 was 3568, checked in by nicolas.vuichard, 8 years ago
modification of the routines where are defined the polygons and their attributes over the grid: single point simulations are now possible, simulations over a single latitudinal or longitudinal band are not. In this later case, a message is printed out
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 38.0 KB

Line
1
2	!! This module define variables for the grid to gathered points.
3	!!
4	!! @call sechiba_main
5	!! @Version : $Revision$, $Date$
6	!!
7	!< $HeadURL$
8	!< $Date$
9	!< $Author$
10	!< $Revision$
11	!!
12	!! @author Marie-Alice Foujols, Jan Polcher and Martial Mancip
13	!!
14	!! This module archives and makes available for all ORCHIDEE routine the information on the grid
15	!! being used. 3 types of grids are foreseen :
16	!! - Regular longitude latitude grid : This is the default and mostly used for global applications.
17	!! - Regular X/Y grid : this is a typical grid for regional models and requires a projection method
18	!! to go from X/y to lon/lat.
19	!! - unstructures grid : This is a general grid where each cell is a polygone. It prepares ORCHIDEE
20	!! for DYNAMICO.
21	!!
22	!! The subroutines have the following role :
23	!! grid_init : this routine will provide the dimensions needed to allocate the memory and the
24	!! characteristics of the grid.
25	!!
26	!! grid_stuff : This subroutine provides the grid details for all land points. Obviously depending
27	!! on the grid type different level of information need to be provided.
28	!!
29	!f90doc MODULEgrid
30	MODULE grid
31
32	USE defprec
33	USE constantes
34	USE mod_orchidee_para
35
36	USE haversine
37	USE module_llxy
38
39	USE ioipsl
40	USE netcdf
41
42	IMPLICIT NONE
43	!
44	!=================================================================================
45	!
46	! Horizontal grid information
47	!
48	!=================================================================================
49	!
50	CHARACTER(LEN=20), SAVE :: GridType
51	!$OMP THREADPRIVATE(GridType)
52	! Describes the grid it can be either of
53	! - RegLonLat
54	! - RegXY
55	! - UnStruct
56	CHARACTER(LEN=20), SAVE :: GridName !! Name of the grid
57	!$OMP THREADPRIVATE(GridName)
58	!
59	INTEGER, SAVE :: NbSegments !! Number of segments in
60	! the polygone defining the grid box
61	!$OMP THREADPRIVATE(NbSegments)
62	!
63	INTEGER, SAVE :: NbNeighb !! Number of neighbours
64	!$OMP THREADPRIVATE(NbNeighb)
65	!
66	! Global map or not.
67	! There is little chance that if iim <=2 and jjm <= 2 that we have global grid.
68	! Furthermore using the second line allows to avoid pole problems for global grids
69	LOGICAL, SAVE :: global = .TRUE.
70	!$OMP THREADPRIVATE(global)
71	!
72	! PARAMETERS
73	! default resolution (m)
74	REAL(r_std), PARAMETER :: default_resolution = 250000.
75	!
76	! VARIABLES
77	!
78	!-
79	!- Variable to help describe the grid
80	!- once the points are gathered.
81	!-
82	!! Limits of the domain
83	REAL(r_std), SAVE :: limit_west, limit_east, &
84	& limit_north, limit_south
85	!$OMP THREADPRIVATE(limit_west, limit_east, limit_north, limit_south)
86	!-
87	!! Geographical coordinates
88	REAL(r_std), ALLOCATABLE, DIMENSION (:,:), SAVE :: lalo
89	!$OMP THREADPRIVATE(lalo)
90	!! index of land points
91	INTEGER, ALLOCATABLE, DIMENSION (:), SAVE :: ilandindex,jlandindex
92	!$OMP THREADPRIVATE(ilandindex, jlandindex)
93	!-
94	!! Fraction of continents.
95	REAL(r_std), ALLOCATABLE, DIMENSION (:), SAVE :: contfrac
96	!$OMP THREADPRIVATE(contfrac)
97	!
98	! indices of the NbNeighb neighbours of each grid point
99	! (1=Northern most vertex and then in clockwise order)
100	! Zero or negative index means that this neighbour is not a land point
101	INTEGER(i_std), ALLOCATABLE, DIMENSION (:,:), SAVE :: neighbours
102	!$OMP THREADPRIVATE(neighbours)
103	!
104	! Heading of the direction out of the grid box either through the vertex
105	! of the mid-segment of the polygon.
106	!
107	REAL(r_std), ALLOCATABLE, DIMENSION(:,:), SAVE :: headings
108	!$OMP THREADPRIVATE(headings)
109	!
110	! Length of segments of the polygon.
111	!
112	REAL(r_std), ALLOCATABLE, DIMENSION(:,:), SAVE :: seglength
113	!$OMP THREADPRIVATE(seglength)
114	!
115	! Area of the grid box
116	!
117	REAL(r_std), ALLOCATABLE, DIMENSION(:), SAVE :: area
118	!$OMP THREADPRIVATE(area)
119	!
120	! Coordinats of the vertices
121	!
122	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: corners
123	!$OMP THREADPRIVATE(corners)
124	!
125	! Resolution remains a temporary variable until the merge of the
126	! re-interfacing of the interpolation by Lluis. One this is done
127	! Resolution will be replaced in the model either by area or seglength.
128	!
129	REAL(r_std), ALLOCATABLE, DIMENSION (:,:), SAVE :: resolution
130	!$OMP THREADPRIVATE(resolution)
131	!
132	!
133	!
134	! Get the direction of the grid
135	!
136	CHARACTER(LEN=2), DIMENSION(2), SAVE, PRIVATE :: grid_dir
137	!$OMP THREADPRIVATE(grid_dir)
138	!
139	INTEGER(i_std), PARAMETER :: MAX_DOMAINS=1
140	!
141	type (proj_info), SAVE, dimension(1:MAX_DOMAINS) :: proj_stack
142	!
143	real(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: dxwrf, dywrf
144	!
145	!
146	!=================================================================================
147	!
148	! Calendar information
149	!
150	!=================================================================================
151	!
152	! The calendar
153	CHARACTER(LEN=20), SAVE :: calendar_str
154	!$OMP THREADPRIVATE(calendar_str)
155	!
156	! The date
157	REAL(r_std), SAVE :: in_julian
158	!$OMP THREADPRIVATE(in_julian)
159	! Diff with day 0
160	REAL(r_std), SAVE :: julian_diff
161	!$OMP THREADPRIVATE(julian_diff)
162	!
163	INTEGER(i_std), SAVE :: year, month, day
164	!$OMP THREADPRIVATE(year, month, day)
165	REAL(r_std), SAVE :: sec
166	!$OMP THREADPRIVATE(sec)
167	!
168	! month_len (d)
169	INTEGER(i_std), SAVE :: month_len
170	!$OMP THREADPRIVATE(month_len)
171	!
172	! year length (d)
173	INTEGER(i_std), SAVE :: year_length=0
174	!$OMP THREADPRIVATE(year_length)
175	!
176	! Ration between calendar year in days (ie 360d or 365d ...) to gregorian year length
177	REAL(r_std), SAVE :: year_spread
178	!$OMP THREADPRIVATE(year_spread)
179	!
180	!
181	INTERFACE grid_tolola
182	MODULE PROCEDURE grid_tolola_scal, grid_tolola_1d, grid_tolola_2d
183	END INTERFACE grid_tolola
184
185	INTERFACE grid_toij
186	MODULE PROCEDURE grid_toij_scal, grid_toij_1d, grid_toij_2d
187	END INTERFACE grid_toij
188	!
189	CONTAINS
190	!
191	!f90doc CONTAINS
192	!
193	!
194	!! =============================================================================================================================
195	!! SUBROUTINE: grid_init
196	!!
197	!>\BRIEF Initialization of grid description distributed by this module to the rest of the model.
198	!!
199	!! DESCRIPTION: Routine which provides the dimension of the grid (number of land points) as well as the
200	!! grid characteristics (type and name) so that the memory can be allocated.
201	!!
202	!! This subroutine is called by intersurf_main_2d or any driver of the model.
203	!!
204	!! \n
205	!_ ==============================================================================================================================
206	!!
207	SUBROUTINE grid_init ( npts, nbseg, gtype, gname, isglobal )
208	!
209	! 0 interface
210	!
211	IMPLICIT NONE
212	!
213	! 0.1 input !
214	!
215	! Domain size
216	INTEGER(i_std), INTENT(in) :: npts !! Number of local continental points
217	INTEGER(i_std), INTENT(in) :: nbseg !! number of segments of the polygone of the mesh
218	CHARACTER(LEN=*), INTENT(in) :: gtype !! Type of grid
219	CHARACTER(LEN=*), INTENT(in) :: gname !! Name of the grid
220	LOGICAL, OPTIONAL :: isglobal
221	!
222	!
223	!
224	CHARACTER(LEN=20) :: gtype_lower
225	!
226	! Verify the information passed and save it in the global variables of the model.
227	!
228	gtype_lower = gtype
229	CALL strlowercase(gtype_lower)
230
231	IF ( INDEX(gtype_lower, "reglonlat") > 0) THEN
232	IF ( nbseg /= 4 ) THEN
233	CALL ipslerr(3, "grid_init", "This regular Lon/lat grid should have 4 segments", &
234	& "per horizontal grid box","")
235	ELSE
236	NbSegments=4
237	ENDIF
238	GridType="RegLonLat"
239	GridName=gridname
240	IF ( PRESENT(isglobal) ) THEN
241	global = isglobal
242	ELSE
243	global = .TRUE.
244	ENDIF
245	ELSE IF ( INDEX(gtype_lower, "regxy") > 0) THEN
246	IF ( nbseg /= 4 ) THEN
247	CALL ipslerr(3, "grid_init", "This regular X/Y grid should have 4 segments", &
248	& "per horizontal grid box","")
249	ELSE
250	NbSegments=4
251	ENDIF
252	GridType="RegXY"
253	GridName=gridname
254	IF ( PRESENT(isglobal) ) THEN
255	global = isglobal
256	ELSE
257	global = .FALSE.
258	ENDIF
259	ELSE IF ( INDEX(gtype_lower, "unstruct") > 0) THEN
260	NbSegments=nbseg
261	GridType="UnStruct"
262	GridName=gridname
263	IF ( PRESENT(isglobal) ) THEN
264	global = isglobal
265	ELSE
266	global = .TRUE.
267	ENDIF
268	ELSE
269	CALL ipslerr(3, "grid_init", "unrecognized grid type.",&
270	& "It has to be either reglatlon, regxy or unstruct","")
271	ENDIF
272	!
273	! Create the internal coordinate table
274	!
275	IF ( (.NOT.ALLOCATED(lalo))) THEN
276	ALLOCATE(lalo(npts,2))
277	lalo(:,:) = val_exp
278	ENDIF
279	!-
280	!- Store variable to help describe the grid
281	!- once the points are gathered.
282	!-
283	NbNeighb=2*NbSegments
284	IF ( (.NOT.ALLOCATED(neighbours))) THEN
285	ALLOCATE(neighbours(npts,NbNeighb))
286	neighbours(:,:) = -999999
287	ENDIF
288	IF ( (.NOT.ALLOCATED(headings))) THEN
289	ALLOCATE(headings(npts,NbNeighb))
290	headings(:,:) = val_exp
291	ENDIF
292	IF ( (.NOT.ALLOCATED(seglength))) THEN
293	ALLOCATE(seglength(npts,NbSegments))
294	seglength(:,:) = val_exp
295	ENDIF
296	IF ( (.NOT.ALLOCATED(corners))) THEN
297	ALLOCATE(corners(npts,NbSegments,2))
298	corners(:,:,:) = val_exp
299	ENDIF
300	IF ( (.NOT.ALLOCATED(area))) THEN
301	ALLOCATE(area(npts))
302	area(:) = val_exp
303	ENDIF
304	!
305	! TEMPORARY
306	!
307	IF ( (.NOT.ALLOCATED(resolution))) THEN
308	ALLOCATE(resolution(npts,2))
309	resolution(:,:) = val_exp
310	ENDIF
311	!
312	!- Store the fraction of the continents only once so that the user
313	!- does not change them afterwards.
314	!
315	IF ( (.NOT.ALLOCATED(contfrac))) THEN
316	ALLOCATE(contfrac(npts))
317	contfrac(:) = val_exp
318	ENDIF
319	!
320	! Allocation of index coordinates ...
321	! JP : these are global fields and should perhaps be allocated somewhere else.
322	IF (.NOT. ALLOCATED(ilandindex)) THEN
323	ALLOCATE(ilandindex(nbp_glo),jlandindex(nbp_glo))
324	ilandindex(:) = -10000000
325	jlandindex(:) = -10000000
326	ENDIF
327	!
328	END SUBROUTINE grid_init
329	!!
330	!!
331	!! =============================================================================================================================
332	!! FUNCTION grid_set
333	!!
334	!>\BRIEF subroutine to set global grid parameters present on all procs
335	!!
336	!! DESCRIPTION:
337	!!
338	!!
339	!!
340	!!
341	!! \n
342	!_ ==============================================================================================================================
343	!!
344	SUBROUTINE grid_set_glo(arg_nbp_lon,arg_nbp_lat,arg_nbp_glo)
345	IMPLICIT NONE
346
347	INTEGER(i_std), INTENT(IN) :: arg_nbp_lon
348	INTEGER(i_std), INTENT(IN) :: arg_nbp_lat
349	INTEGER(i_std), INTENT(IN),OPTIONAL :: arg_nbp_glo
350	iim_g=arg_nbp_lon
351	jjm_g=arg_nbp_lat
352	IF (PRESENT(arg_nbp_glo)) nbp_glo=arg_nbp_glo
353	END SUBROUTINE grid_set_glo
354	!! =============================================================================================================================
355	!! FUNCTION grid_set/allocate_glo
356	!!
357	!>\BRIEF subroutines to allocate variables present on all procs
358	!!
359	!! DESCRIPTION:
360	!!
361	!!
362	!!
363	!!
364	!! \n
365	!_ ==============================================================================================================================
366	!!
367	SUBROUTINE grid_allocate_glo(nbseg)
368	!
369	IMPLICIT NONE
370	! 0.1 input !
371	!
372	! Domain size
373	INTEGER(i_std), INTENT(in) :: nbseg !! number of segments of the polygone of the mesh
374	!
375	! In case the allocation of the grid is called before the initialisation,
376	! we already set the number of segments.
377	! This will be done properly in grid_init.
378	!
379	IF ( NbSegments < 3 ) THEN
380	NbSegments = nbseg
381	NbNeighb=2*NbSegments
382	ENDIF
383	!
384	!
385	ALLOCATE(neighbours_g(nbp_glo,NbNeighb))
386	ALLOCATE(headings_g(nbp_glo,NbNeighb))
387	ALLOCATE(seglength_g(nbp_glo,NbSegments))
388	ALLOCATE(corners_g(nbp_glo,NbSegments,2))
389	ALLOCATE(area_g(nbp_glo))
390	!
391	! TEMPORARY
392	!
393	ALLOCATE(resolution_g(nbp_glo,2))
394	!
395	! Allocate other variables
396	!
397	ALLOCATE(lalo_g(nbp_glo,2), contfrac_g(nbp_glo),index_g(nbp_glo))
398	ALLOCATE(lon_g(iim_g, jjm_g), lat_g(iim_g, jjm_g), zlev_g(iim_g, jjm_g))
399	!
400	END SUBROUTINE grid_allocate_glo
401	!!
402	!! =============================================================================================================================
403	!! SUBROUTINE: grid_stuff
404	!!
405	!>\BRIEF transfers the global horizontal grid information to ORCHIDEE in the case of grid regular in Longitude
406	!! and Latitude.
407	!!
408	!! DESCRIPTION:
409	!!
410	!!
411	!! This subroutine is called by intersurf_main_2d or any driver of the model.
412	!!
413	!! \n
414	!_ ==============================================================================================================================
415	!!
416	SUBROUTINE grid_stuff (npts_glo, iim, jjm, grid_lon, grid_lat, kindex, contfrac_tmp)
417	!
418	! 0 interface
419	!
420	IMPLICIT NONE
421	!
422	! 0.1 input !
423
424	! Domain size
425	INTEGER(i_std), INTENT(in) :: npts_glo
426	! Size of cartesian grid
427	INTEGER(i_std), INTENT(in) :: iim, jjm
428	! Longitudes on cartesian grid
429	REAL(r_std), DIMENSION(iim,jjm), INTENT(in) :: grid_lon
430	! Latitudes on cartesian grid
431	REAL(r_std), DIMENSION(iim,jjm), INTENT(in) :: grid_lat
432	! Index of land point on 2D map (in local position)
433	INTEGER(i_std), DIMENSION(npts_glo), INTENT(in) :: kindex
434	! The fraction of continent in the grid box [0-1]
435	REAL(r_std), DIMENSION(npts_glo), OPTIONAL, INTENT(in) :: contfrac_tmp
436	!
437	!
438	! =========================================================================
439
440	IF ( printlev >= 4 ) WRITE(numout,*) 'Entering grid_stuff'
441
442	! default resolution
443	IF ( printlev >=2 ) WRITE(numout,*) 'grid stuff: default resolution (m): ',default_resolution
444	!
445	!-
446	IF (is_root_prc) THEN
447	!
448	CALL grid_topolylist(GridType, NbSegments, npts_glo, iim, jjm, grid_lon, grid_lat, kindex, &
449	& global, corners_g, neighbours_g, headings_g, seglength_g, area_g, ilandindex, jlandindex)
450	!
451	IF (PRESENT(contfrac_tmp)) THEN
452	!
453	! Transfer the contfrac into the array managed in this module.
454	!
455	contfrac_g(:) = contfrac_tmp(:)
456	ENDIF
457	!
458	ENDIF
459	!
460	! With this the description of the grid is complete and the information
461	! can be scattered to all processors.
462	!
463	CALL grid_scatter()
464	!
465	IF ( printlev >= 3 ) WRITE(numout,*) 'Leaving grid_stuff'
466
467	END SUBROUTINE grid_stuff
468	!!
469	!! =============================================================================================================================
470	!! SUBROUTINE: grid_topolylist
471	!!
472	!>\BRIEF This routine transforms a regular grid into a list of polygons which are defined by the following
473	!! quantities :
474	!!
475	!! corners : the n vertices of the polugon in longitude and latitude
476	!! neighbours : the neighbouring land grid box for each of the vertices and segments
477	!! headings : the direction in which the neighbour is
478	!! seglength : the lenght of each segment
479	!! area : the area of the polygon
480	!! ilindex, jlindex : provides the i,j coordinates of the mesh in the global grid.
481	!!
482	!! DESCRIPTION:
483	!!
484	!! \n
485	!_ ==============================================================================================================================
486	!!
487	SUBROUTINE grid_topolylist(gtype, nbseg, nland, iim, jjm, grid_lon, grid_lat, kindex, &
488	& globalg, corners_loc, neighbours_loc, headings_loc, seglength_loc, &
489	& area_loc, ilindex_loc, jlindex_loc)
490	!
491	! 0 interface
492	!
493	IMPLICIT NONE
494	!
495	! 0.1 input !
496	! Grid type
497	CHARACTER(LEN=20), INTENT(in) :: gtype
498	! Number of segments for each polygon
499	INTEGER(i_std), INTENT(in) :: nbseg
500	! Number of land points on the grid
501	INTEGER(i_std), INTENT(in) :: nland
502	! Size of cartesian grid
503	INTEGER(i_std), INTENT(in) :: iim, jjm
504	! Longitudes on cartesian grid
505	REAL(r_std), DIMENSION(iim,jjm), INTENT(in) :: grid_lon
506	! Latitudes on cartesian grid
507	REAL(r_std), DIMENSION(iim,jjm), INTENT(in) :: grid_lat
508	! Index of land point on 2D map (in local position)
509	INTEGER(i_std), DIMENSION(nland), INTENT(in) :: kindex
510	!
511	! 0.2 Output
512	!
513	LOGICAL, INTENT(out) :: globalg
514	!
515	REAL(r_std), DIMENSION(nland,nbseg,2), INTENT(out) :: corners_loc
516	INTEGER(i_std), DIMENSION(nland,nbseg*2), INTENT(out) :: neighbours_loc
517	REAL(r_std), DIMENSION(nland,nbseg*2), INTENT(out) :: headings_loc
518	REAL(r_std), DIMENSION(nland,nbseg), INTENT(out) :: seglength_loc
519	REAL(r_std), DIMENSION(nland), INTENT(out) :: area_loc
520	INTEGER(i_std), DIMENSION(nland), INTENT(out) :: ilindex_loc, jlindex_loc
521	!
522	! 0.3 Local variables
523	!
524	INTEGER(i_std) :: i, is, iss
525	REAL(r_std), DIMENSION(nland,2) :: center
526	REAL(r_std) :: maxdellon, mindellon, maxlon, minlon
527	REAL(r_std), DIMENSION(nland,nbseg*2) :: lonpoly, latpoly
528	!
529	IF ( INDEX(gtype,"RegLonLat") > 0 ) THEN
530	!
531	! If we are in regular Lon Lat, then we test just the longitude and see if we span 0-360deg.
532	!
533	maxdellon=MAXVAL(ABS(grid_lon(1:iim-1,1)-grid_lon(2:iim,1)))
534	mindellon=MINVAL(ABS(grid_lon(1:iim-1,1)-grid_lon(2:iim,1)))
535	maxlon=MAXVAL(grid_lon(1:iim,1))
536	minlon=MINVAL(grid_lon(1:iim,1))
537	!
538	! test if it could be a global grid on 0 -> 360
539	!
540	IF ( minlon > 0 .AND. maxlon > 180 ) THEN
541	IF ( (minlon - maxdellon/2.0 ) <= 0 .AND. (maxlon + maxdellon/2.0) >= 360) THEN
542	globalg = .TRUE.
543	ELSE
544	globalg = .FALSE.
545	ENDIF
546	!
547	! Test if it could be a -180 to 180 grid
548	!
549	ELSE IF ( minlon < 0 .AND. maxlon > 0 ) THEN
550	IF ( (minlon - maxdellon/2.0 ) <= -180 .AND. (maxlon + maxdellon/2.0) >= 180) THEN
551	globalg = .TRUE.
552	ELSE
553	globalg = .FALSE.
554	ENDIF
555	!
556	! If neither condition is met then it cannot be global.
557	!
558	ELSE
559	globalg = .FALSE.
560	ENDIF
561	ELSE IF ( gtype == "RegXY" ) THEN
562	!
563	! The hypothesis is that if we are in RegXY then we are not global
564	!
565	globalg = .FALSE.
566	ELSE
567	STOP "Unknown grid"
568	ENDIF
569	!
570	! 2.0 Transform the grid into a list of polygones while keeping the neighbour relations
571	! between these polygones.
572	!
573	! Each polygone starts with a vertex and alternates vertices and mid-points of segments.
574	!
575	IF (nland == 1) THEN
576	CALL haversine_singlepointploy(iim, jjm, grid_lon, grid_lat, nland, kindex, global, &
577	& nbseg, lonpoly, latpoly, center, &
578	& neighbours_loc, ilindex_loc, jlindex_loc)
579	ELSE IF ( INDEX(gtype, "RegLonLat") > 0 ) THEN
580	CALL haversine_reglatlontoploy(iim, jjm, grid_lon, grid_lat, nland, kindex, global, &
581	& nbseg, lonpoly, latpoly, center, &
582	& neighbours_loc, ilindex_loc, jlindex_loc)
583	ELSE IF ( INDEX(gtype, "RegXY") > 0 ) THEN
584	CALL haversine_regxytoploy(iim, jjm, grid_lon, grid_lat, nland, kindex, proj_stack, &
585	& nbseg, lonpoly, latpoly, center, &
586	& neighbours_loc, ilindex_loc, jlindex_loc)
587	ELSE
588	STOP "Unknown grid"
589	ENDIF
590	!
591	! Save the longitude and latitudes nbseg corners (=vertices) of the polygones
592	!
593	DO i=1,nland
594	DO is=1,nbseg
595	iss=(is-1)*2+1
596	corners_loc(i,is,1) = lonpoly(i,iss)
597	corners_loc(i,is,2) = latpoly(i,iss)
598	ENDDO
599	ENDDO
600	!
601	! Get the heading normal to the 4 segments and through the 4 corners.
602	!
603	CALL haversine_polyheadings(nland, nbseg, lonpoly, latpoly, center, headings_loc)
604	!
605	! Order the points of the polygone in clockwise order Starting with the northern most
606	!
607	CALL haversine_polysort(nland, nbseg, lonpoly, latpoly, headings_loc, neighbours_loc)
608	!
609	! Compute the segment length and area.
610	! For the RegLonLat we have specific calculations for seglength and area.
611	! For projected regular grids we use the great cicle assumption for the segments
612	! but the projected area.
613	! For unstructured grid we use the most general routines.
614	!
615	IF ( INDEX(gtype, "RegLonLat") > 0 ) THEN
616	CALL haversine_laloseglen(nland, nbseg, lonpoly, latpoly, seglength_loc)
617	CALL haversine_laloarea(nland, nbseg, seglength_loc, area_loc)
618	ELSE IF ( INDEX(gtype, "RegXY") > 0 ) THEN
619	CALL haversine_polyseglen(nland, nbseg, lonpoly, latpoly, seglength_loc)
620	CALL haversine_xyarea(nland, nbseg, ilindex_loc, jlindex_loc, dxwrf, dywrf, area_loc)
621	ELSE
622	CALL haversine_polyseglen(nland, nbseg, lonpoly, latpoly, seglength_loc)
623	CALL haversine_polyarea(nland, nbseg, lonpoly, latpoly, area_loc)
624	ENDIF
625	! Compute the area
626
627	!
628	END SUBROUTINE grid_topolylist
629	!!
630	!!
631	!!
632	!! =============================================================================================================================
633	!! SUBROUTINE: grid_scatter
634	!!
635	!>\BRIEF Scatter the grid information so that each processor knows the characteristics of the grid it works on.
636	!!
637	!! DESCRIPTION:
638	!!
639	!!
640	!! The grid information has been computed for the entire grid on the root processor. Now we give each processor
641	!! the information of the piece of the grid it works on. This concerns the following variables describing the grid :
642	!! - area
643	!! - resolution
644	!! - neighbours
645	!! - contfrac : fraction of continent
646	!!
647	!! Should ilandindex and jlandindex not b initialized, we catch-up here. This field is the same on all processors.
648	!!
649	!! TODO :
650	!! This code should get the grid describing fields as arguments and then writem into the *_g variables on
651	!! root_prc before scattering. This would allow to compute the grid characteristics in any subroutine
652	!! fore calling grid_scatter.
653	!!
654	!!
655	!!
656	!! \n
657	!_ ==============================================================================================================================
658	!!
659	!!
660	SUBROUTINE grid_scatter()
661	!
662	!
663	INTEGER(i_std) :: i, ip, jp
664	!
665	IF ( MAXVAL(ilandindex) < 0 .AND. MAXVAL(jlandindex) < 0 ) THEN
666	DO i = 1, nbp_glo
667	!
668	! 1 find numbers of the latitude and longitude of each point
669	!
670
671	! index of latitude
672	jp = INT( (index_g(i)-1) /iim_g ) + 1
673
674	! index of longitude
675	ip = index_g(i) - ( jp-1 ) * iim_g
676	!
677	! Save this information for usage in other modules.
678	!
679	ilandindex(i)=ip
680	jlandindex(i)=jp
681	!
682	ENDDO
683	ENDIF
684	!
685	CALL scatter(neighbours_g, neighbours)
686	CALL scatter(contfrac_g, contfrac)
687	CALL scatter(headings_g, headings)
688	CALL scatter(seglength_g, seglength)
689	CALL scatter(corners_g, corners)
690	CALL scatter(area_g, area)
691	!
692	! TEMPORARY section for resolution
693	!
694	IF ( is_root_prc) THEN
695	IF ( INDEX(GridType,"Reg") > 0 ) THEN
696	resolution_g(:,1) = (seglength_g(:,1)+seglength_g(:,3))/2.0
697	resolution_g(:,2) = (seglength_g(:,2)+seglength_g(:,4))/2.0
698	ELSE
699	CALL ipslerr(3, "grid_scatter", "unsupported grid type.",&
700	& "As long as resolution has not been replaced,",&
701	& "ORCHIDEE cannot run on anything other than regular grids.")
702	ENDIF
703	ENDIF
704	CALL scatter(resolution_g, resolution)
705	!
706	!
707	IF ( printlev >=4 ) THEN
708	WRITE(numout,*) 'grid_scatter > seglength = ', seglength(1,:)
709	WRITE(numout,*) 'grid_scatter > neighbours = ', neighbours(1,:)
710	WRITE(numout,*) 'grid_scatter > contfrac = ', contfrac(1)
711	WRITE(numout,*) 'grid_scatter > area = ', area(1)
712	ENDIF
713	!
714	END SUBROUTINE grid_scatter
715	!!
716	!!
717	!! =============================================================================================================================
718	!! SUBROUTINE: grid_initproj
719	!!
720	!>\BRIEF Routine to initialise the projection
721	!!
722	!! DESCRIPTION:
723	!!
724	!!
725	!! This subroutine is called by the routine whichs ets-up th grid on which ORCHIDEE is to run.
726	!! The aim is to set-upu the projection so that all the grid variables needed by ORCHIDEE can
727	!! be computed in grid_stuff_regxy
728	!!
729	!! \n
730	!_ ==============================================================================================================================
731	!!
732	!!
733	SUBROUTINE grid_initproj (fid, iim, jjm)
734	!
735	!
736	! 0 interface
737	!
738	IMPLICIT NONE
739	!
740	! 0.1 input !
741	!
742	! Domain size
743	INTEGER(i_std), INTENT(in) :: fid
744	INTEGER(i_std), INTENT(in) :: iim, jjm
745	!
746	! 0.2 Local variables
747	!
748	INTEGER(i_std) :: current_proj, idom, iret, lonid, latid, numLons, numLats
749	INTEGER, DIMENSION(nf90_max_var_dims) :: dimIDs
750	REAL(r_std) :: user_stand_lon, user_truelat1, user_truelat2, user_dxkm, user_dykm
751	REAL(r_std) :: user_dlat, user_dlon, user_known_x, user_known_y, user_known_lat, user_known_lon
752	REAL(r_std), DIMENSION(16) :: corner_lons, corner_lats
753	!
754	INTEGER(i_std) :: iv, i, j
755	CHARACTER(LEN=20) :: varname
756	REAL(r_std) :: dx, dy, dtx, dty, coslat
757	REAL(r_std), ALLOCATABLE, DIMENSION (:) :: LON, LAT
758	REAL(r_std), ALLOCATABLE, DIMENSION (:,:) :: mapfac_x, mapfac_y
759	!
760	!
761	! Only one domain is possible for the moment
762	!
763	idom=1
764	CALL map_init(proj_stack(idom))
765	!
766	! Does ORCHIDEE have the same Earth Radius as the map projection ?
767	!
768	IF ( ABS(R_Earth-EARTH_RADIUS_M) > 0.1 ) THEN
769	WRITE(,) "Earth Radius in WRF : ", EARTH_RADIUS_M
770	WRITE(,) "Earth Radius in ORCHIDEE : ", R_Earth
771	CALL ipslerr (3,'grid_initproj','The Earth radius is not the same in the projection module and ORCHIDEE',&
772	& " ", " ")
773	ENDIF
774	!
775	! Get parameters of the projection from the netCDF file
776	!
777	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "MAP_PROJ", current_proj)
778	!
779	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "STAND_LON", user_stand_lon)
780	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "TRUELAT1", user_truelat1)
781	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "TRUELAT2", user_truelat2)
782	!
783	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "DX", user_dxkm)
784	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "DY", user_dykm)
785	user_dlat = undef
786	user_dlon = undef
787	!
788	IF ( current_proj == PROJ_LATLON ) THEN
789	!
790	iret = NF90_inq_VARID(fid, "XLONG_M",lonid)
791	iret = NF90_INQUIRE_VARIABLE(fid, lonid, dimids = dimIDs)
792	iret = NF90_INQUIRE_DIMENSION(fid, dimIDs(1), len = numLons)
793	iret = NF90_INQUIRE_DIMENSION(fid, dimIDs(2), len = numLats)
794	ALLOCATE(LON(numLons))
795	iret = NF90_GET_VAR(fid, lonid, LON(:), start = (/ 1, 1, 1 /), count = (/ numLons, 1, 1 /))
796
797	iret = NF90_inq_VARID(fid, "XLAT_M",latid)
798	ALLOCATE(LAT(numLats))
799	iret = NF90_GET_VAR(fid, latid, LAT(:), start = (/ 1, 1, 1 /), count = (/ 1, numLats, 1 /))
800
801	user_dlon = (LON(numLons) - LON(1)) / (numLons - 1)
802	user_dlat = (LAT(numLats) - LAT(1)) / (numLats - 1)
803
804	DEALLOCATE(LON,LAT)
805
806	ENDIF
807	! Unable to know from where to get the information
808	user_known_x = 1
809	user_known_y = 1
810	!
811	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "corner_lats", corner_lats)
812	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "corner_lons", corner_lons)
813	user_known_lat = corner_lats(1)
814	user_known_lon = corner_lons(1)
815	!
816	! Read mapfactor, land mask and orography
817	!
818	!
819	! Allocation
820	!
821	ALLOCATE(mapfac_x(iim,jjm))
822	ALLOCATE(mapfac_y(iim,jjm))
823	ALLOCATE(dxwrf(iim,jjm))
824	ALLOCATE(dywrf(iim,jjm))
825	!
826	varname = "MAPFAC_MX"
827	iret = NF90_INQ_VARID (fid, varname, iv)
828	IF (iret /= NF90_NOERR) THEN
829	CALL ipslerr (3,'WRFdomain_Read',"Could not find variable ", varname," ")
830	ELSE
831	iret = NF90_GET_VAR (fid,iv,mapfac_x)
832	ENDIF
833	varname = "MAPFAC_MY"
834	iret = NF90_INQ_VARID (fid, varname, iv)
835	IF (iret /= NF90_NOERR) THEN
836	CALL ipslerr (3,'WRFdomain_Read',"Could not find variable ", varname," ")
837	ELSE
838	iret = NF90_GET_VAR (fid,iv,mapfac_y)
839	ENDIF
840	!
841	! Initilize the projection
842	!
843	if (current_proj == PROJ_LATLON) then
844	call map_set(current_proj, proj_stack(idom), &
845	lat1=user_known_lat, &
846	lon1=user_known_lon, &
847	knowni=user_known_x, &
848	knownj=user_known_y, &
849	latinc=user_dlat, &
850	loninc=user_dlon, &
851	r_earth=R_Earth)
852
853	else if (current_proj == PROJ_MERC) then
854	call map_set(current_proj, proj_stack(idom), &
855	truelat1=user_truelat1, &
856	lat1=user_known_lat, &
857	lon1=user_known_lon, &
858	knowni=user_known_x, &
859	knownj=user_known_y, &
860	dx=user_dxkm, &
861	r_earth=R_Earth)
862
863	else if (current_proj == PROJ_CYL) then
864	call ipslerr(3,"grid_initproj",'Should not have PROJ_CYL as projection for',&
865	'source data in push_source_projection()', " ")
866
867	else if (current_proj == PROJ_CASSINI) then
868	call ipslerr(3,"grid_initproj",'Should not have PROJ_CASSINI as projection for', &
869	'source data in push_source_projection()', " ")
870
871	else if (current_proj == PROJ_LC) then
872	call map_set(current_proj, proj_stack(idom), &
873	truelat1=user_truelat1, &
874	truelat2=user_truelat2, &
875	stdlon=user_stand_lon, &
876	lat1=user_known_lat, &
877	lon1=user_known_lon, &
878	knowni=user_known_x, &
879	knownj=user_known_y, &
880	dx=user_dxkm, &
881	r_earth=R_Earth)
882
883	else if (current_proj == PROJ_ALBERS_NAD83) then
884	call map_set(current_proj, proj_stack(idom), &
885	truelat1=user_truelat1, &
886	truelat2=user_truelat2, &
887	stdlon=user_stand_lon, &
888	lat1=user_known_lat, &
889	lon1=user_known_lon, &
890	knowni=user_known_x, &
891	knownj=user_known_y, &
892	dx=user_dxkm, &
893	r_earth=R_Earth)
894
895	else if (current_proj == PROJ_PS) then
896	call map_set(current_proj, proj_stack(idom), &
897	truelat1=user_truelat1, &
898	stdlon=user_stand_lon, &
899	lat1=user_known_lat, &
900	lon1=user_known_lon, &
901	knowni=user_known_x, &
902	knownj=user_known_y, &
903	dx=user_dxkm, &
904	r_earth=R_Earth)
905
906	else if (current_proj == PROJ_PS_WGS84) then
907	call map_set(current_proj, proj_stack(idom), &
908	truelat1=user_truelat1, &
909	stdlon=user_stand_lon, &
910	lat1=user_known_lat, &
911	lon1=user_known_lon, &
912	knowni=user_known_x, &
913	knownj=user_known_y, &
914	dx=user_dxkm, &
915	r_earth=R_Earth)
916
917	else if (current_proj == PROJ_GAUSS) then
918	call map_set(current_proj, proj_stack(idom), &
919	lat1=user_known_lat, &
920	lon1=user_known_lon, &
921	nlat=nint(user_dlat), &
922	loninc=user_dlon, &
923	r_earth=R_Earth)
924
925	else if (current_proj == PROJ_ROTLL) then
926	call ipslerr(3 ,"grid_initproj",'Should not have PROJ_ROTLL as projection for', &
927	'source data in push_source_projection() as not yet implemented', '')
928	end if
929	!
930	! Transform the mapfactors into dx and dy to be used for the description of the polygons and
931	! interpolations.
932	!
933	DO i=1,iim
934	DO j=1,jjm
935	!
936	IF (proj_stack(idom)%code /= PROJ_LATLON ) THEN
937	dx = proj_stack(idom)%dx
938	! Some projections in WRF do not store dy, in that case dy=dx.
939	IF ( proj_stack(idom)%dy > 0 ) THEN
940	dy = proj_stack(idom)%dy
941	ELSE
942	dy = proj_stack(idom)%dx
943	ENDIF
944	dxwrf(i,j) = dx/mapfac_x(i,j)
945	dywrf(i,j) = dy/mapfac_y(i,j)
946	ELSE
947	!
948	! The LatLon projection is also a special case as here it is not the dx and dy
949	! which are stored in the projection file but the increments in Lon and Lat.
950	!
951	dtx = proj_stack(idom)%loninc
952	dty = proj_stack(idom)%latinc
953	coslat = COS(lat(j) * pi/180. )
954	dxwrf(i,j) = dtx * pi/180. * R_Earth * coslat
955	dywrf(i,j) = dty * pi/180. * R_Earth
956	!
957	ENDIF
958	!
959	ENDDO
960	ENDDO
961	!
962	END SUBROUTINE grid_initproj
963	!
964	!
965	!
966	!=========================================================================================
967	!
968	SUBROUTINE grid_tolola_scal (ri, rj, lon, lat)
969	!
970	!
971	! Argument
972	REAL(r_std), INTENT(in) :: ri, rj
973	REAL(r_std), INTENT(out) :: lon, lat
974	!
975	!
976	IF ( proj_stack(1)%code < undef_int ) THEN
977	!
978	CALL ij_to_latlon(proj_stack(1), ri, rj, lat, lon)
979	!
980	ELSE
981	CALL ipslerr(3, "grid_tolola_scal", "Projection not initilized"," "," ")
982	ENDIF
983	!
984	END SUBROUTINE grid_tolola_scal
985	!
986	!=========================================================================================
987	!
988	SUBROUTINE grid_tolola_1d (ri, rj, lon, lat)
989	!
990	!
991	! Argument
992	REAL(r_std), INTENT(in), DIMENSION(:) :: ri, rj
993	REAL(r_std), INTENT(out), DIMENSION(:) :: lon, lat
994	!
995	! Local
996	INTEGER :: i, imax
997	!
998	imax=SIZE(lon)
999	!
1000	IF ( proj_stack(1)%code < undef_int ) THEN
1001	DO i=1,imax
1002	!
1003	CALL ij_to_latlon(proj_stack(1), ri(i), rj(i), lat(i), lon(i))
1004	!
1005	ENDDO
1006	ELSE
1007	CALL ipslerr(3, "grid_tolola_1d", "Projection not initilized"," "," ")
1008	ENDIF
1009	!
1010	END SUBROUTINE grid_tolola_1d
1011	!
1012	!=========================================================================================
1013	!
1014	SUBROUTINE grid_tolola_2d (ri, rj, lon, lat)
1015	!
1016	!
1017	! Argument
1018	REAL(r_std), INTENT(in), DIMENSION(:,:) :: ri, rj
1019	REAL(r_std), INTENT(out), DIMENSION(:,:) :: lon, lat
1020	!
1021	! Local
1022	INTEGER :: i, imax, j, jmax
1023	!
1024	imax=SIZE(lon,DIM=1)
1025	jmax=SIZE(lon,DIM=2)
1026	!
1027	IF ( proj_stack(1)%code < undef_int ) THEN
1028	DO i=1,imax
1029	DO j=1,jmax
1030	!
1031	CALL ij_to_latlon(proj_stack(1), ri(i,j), rj(i,j), lat(i,j), lon(i,j))
1032	!
1033	ENDDO
1034	ENDDO
1035	ELSE
1036	CALL ipslerr(3, "grid_tolola_2d", "Projection not initilized"," "," ")
1037	ENDIF
1038	!
1039	END SUBROUTINE grid_tolola_2d
1040	!
1041	!=========================================================================================
1042	!
1043	SUBROUTINE grid_toij_scal (lon, lat, ri, rj)
1044	!
1045	!
1046	! Argument
1047	REAL(r_std), INTENT(in) :: lon, lat
1048	REAL(r_std), INTENT(out) :: ri, rj
1049	!
1050	!
1051	IF ( proj_stack(1)%code < undef_int ) THEN
1052	!
1053	CALL latlon_to_ij(proj_stack(1), lat, lon, ri, rj)
1054	!
1055	ELSE
1056	CALL ipslerr(3, "grid_toij_scal", "Projection not initilized"," "," ")
1057	ENDIF
1058	!
1059	END SUBROUTINE grid_toij_scal
1060	!
1061	!=========================================================================================
1062	!
1063	SUBROUTINE grid_toij_1d (lon, lat, ri, rj)
1064	!
1065	!
1066	! Argument
1067	REAL(r_std), INTENT(in), DIMENSION(:) :: lon, lat
1068	REAL(r_std), INTENT(out), DIMENSION(:) :: ri, rj
1069	!
1070	! Local
1071	INTEGER :: i, imax
1072	!
1073	imax=SIZE(lon)
1074	!
1075	IF ( proj_stack(1)%code < undef_int ) THEN
1076	DO i=1,imax
1077	!
1078	CALL latlon_to_ij(proj_stack(1), lat(i), lon(i), ri(i), rj(i))
1079	!
1080	ENDDO
1081	ELSE
1082	CALL ipslerr(3, "grid_toij_1d", "Projection not initilized"," "," ")
1083	ENDIF
1084	!
1085	END SUBROUTINE grid_toij_1d
1086	!
1087	!=========================================================================================
1088	!
1089	SUBROUTINE grid_toij_2d (lon, lat, ri, rj)
1090	!
1091	!
1092	! Argument
1093	REAL(r_std), INTENT(in), DIMENSION(:,:) :: lon, lat
1094	REAL(r_std), INTENT(out), DIMENSION(:,:) :: ri, rj
1095	!
1096	! Local
1097	INTEGER :: i, imax, j, jmax
1098	!
1099	imax=SIZE(lon,DIM=1)
1100	jmax=SIZE(lon,DIM=2)
1101	!
1102	IF ( proj_stack(1)%code < undef_int ) THEN
1103	DO i=1,imax
1104	DO j=1,jmax
1105	!
1106	CALL latlon_to_ij(proj_stack(1), lat(i,j), lon(i,j), ri(i,j), rj(i,j))
1107	!
1108	ENDDO
1109	ENDDO
1110	ELSE
1111	CALL ipslerr(3, "grid_toij_2d", "Projection not initilized"," "," ")
1112	ENDIF
1113	!
1114	END SUBROUTINE grid_toij_2d
1115	!
1116	!
1117	!=========================================================================================
1118	!
1119	!
1120	END MODULE grid

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

Download in other formats: