Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-57/src_global/grid.f90 @ 5143

Last change on this file since 5143 was 3965, checked in by jan.polcher, 8 years ago

Merge with trunk at revision3959.
This includes all the developments made for CMIP6 and passage to XIOS2.
All conflicts are resolved and the code compiles.

But it still does not link because of an "undefined reference to `_intel_fast_memmove'"

Property svn:keywords set to Revision Date HeadURL Date Author Revision

File size: 37.9 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	USE grid_var
43
44	IMPLICIT NONE
45	!
46	!=================================================================================
47	!
48	! Horizontal grid information
49	!
50	!=================================================================================
51	!
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	INTEGER(i_std) :: i, ip, jp
663	!
664	IF ( MAXVAL(ilandindex) < 0 .AND. MAXVAL(jlandindex) < 0 ) THEN
665	DO i = 1, nbp_glo
666	!
667	! 1 find numbers of the latitude and longitude of each point
668	!
669
670	! index of latitude
671	jp = INT( (index_g(i)-1) /iim_g ) + 1
672
673	! index of longitude
674	ip = index_g(i) - ( jp-1 ) * iim_g
675	!
676	! Save this information for usage in other modules.
677	!
678	ilandindex(i)=ip
679	jlandindex(i)=jp
680	!
681	ENDDO
682	ENDIF
683	!
684	CALL scatter(neighbours_g, neighbours)
685	CALL scatter(contfrac_g, contfrac)
686	CALL scatter(headings_g, headings)
687	CALL scatter(seglength_g, seglength)
688	CALL scatter(corners_g, corners)
689	CALL scatter(area_g, area)
690	!
691	! TEMPORARY section for resolution
692	!
693	IF ( is_root_prc) THEN
694	IF ( INDEX(GridType,"Reg") > 0 ) THEN
695	resolution_g(:,1) = (seglength_g(:,1)+seglength_g(:,3))/2.0
696	resolution_g(:,2) = (seglength_g(:,2)+seglength_g(:,4))/2.0
697	ELSE
698	CALL ipslerr(3, "grid_scatter", "unsupported grid type.",&
699	& "As long as resolution has not been replaced,",&
700	& "ORCHIDEE cannot run on anything other than regular grids.")
701	ENDIF
702	ENDIF
703	CALL scatter(resolution_g, resolution)
704	!
705	!
706	IF ( printlev >=4 ) THEN
707	WRITE(numout,*) 'grid_scatter > seglength = ', seglength(1,:)
708	WRITE(numout,*) 'grid_scatter > neighbours = ', neighbours(1,:)
709	WRITE(numout,*) 'grid_scatter > contfrac = ', contfrac(1)
710	WRITE(numout,*) 'grid_scatter > area = ', area(1)
711	ENDIF
712	!
713	END SUBROUTINE grid_scatter
714	!!
715	!!
716	!! =============================================================================================================================
717	!! SUBROUTINE: grid_initproj
718	!!
719	!>\BRIEF Routine to initialise the projection
720	!!
721	!! DESCRIPTION:
722	!!
723	!!
724	!! This subroutine is called by the routine whichs ets-up th grid on which ORCHIDEE is to run.
725	!! The aim is to set-upu the projection so that all the grid variables needed by ORCHIDEE can
726	!! be computed in grid_stuff_regxy
727	!!
728	!! \n
729	!_ ==============================================================================================================================
730	!!
731	!!
732	SUBROUTINE grid_initproj (fid, iim, jjm)
733	!
734	!
735	! 0 interface
736	!
737	IMPLICIT NONE
738	!
739	! 0.1 input !
740	!
741	! Domain size
742	INTEGER(i_std), INTENT(in) :: fid
743	INTEGER(i_std), INTENT(in) :: iim, jjm
744	!
745	! 0.2 Local variables
746	!
747	INTEGER(i_std) :: current_proj, idom, iret, lonid, latid, numLons, numLats
748	INTEGER, DIMENSION(nf90_max_var_dims) :: dimIDs
749	REAL(r_std) :: user_stand_lon, user_truelat1, user_truelat2, user_dxkm, user_dykm
750	REAL(r_std) :: user_dlat, user_dlon, user_known_x, user_known_y, user_known_lat, user_known_lon
751	REAL(r_std), DIMENSION(16) :: corner_lons, corner_lats
752	!
753	INTEGER(i_std) :: iv, i, j
754	CHARACTER(LEN=20) :: varname
755	REAL(r_std) :: dx, dy, dtx, dty, coslat
756	REAL(r_std), ALLOCATABLE, DIMENSION (:) :: LON, LAT
757	REAL(r_std), ALLOCATABLE, DIMENSION (:,:) :: mapfac_x, mapfac_y
758	!
759	!
760	! Only one domain is possible for the moment
761	!
762	idom=1
763	CALL map_init(proj_stack(idom))
764	!
765	! Does ORCHIDEE have the same Earth Radius as the map projection ?
766	!
767	IF ( ABS(R_Earth-EARTH_RADIUS_M) > 0.1 ) THEN
768	WRITE(,) "Earth Radius in WRF : ", EARTH_RADIUS_M
769	WRITE(,) "Earth Radius in ORCHIDEE : ", R_Earth
770	CALL ipslerr (3,'grid_initproj','The Earth radius is not the same in the projection module and ORCHIDEE',&
771	& " ", " ")
772	ENDIF
773	!
774	! Get parameters of the projection from the netCDF file
775	!
776	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "MAP_PROJ", current_proj)
777	!
778	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "STAND_LON", user_stand_lon)
779	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "TRUELAT1", user_truelat1)
780	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "TRUELAT2", user_truelat2)
781	!
782	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "DX", user_dxkm)
783	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "DY", user_dykm)
784	user_dlat = undef
785	user_dlon = undef
786	!
787	IF ( current_proj == PROJ_LATLON ) THEN
788	!
789	iret = NF90_inq_VARID(fid, "XLONG_M",lonid)
790	iret = NF90_INQUIRE_VARIABLE(fid, lonid, dimids = dimIDs)
791	iret = NF90_INQUIRE_DIMENSION(fid, dimIDs(1), len = numLons)
792	iret = NF90_INQUIRE_DIMENSION(fid, dimIDs(2), len = numLats)
793	ALLOCATE(LON(numLons))
794	iret = NF90_GET_VAR(fid, lonid, LON(:), start = (/ 1, 1, 1 /), count = (/ numLons, 1, 1 /))
795
796	iret = NF90_inq_VARID(fid, "XLAT_M",latid)
797	ALLOCATE(LAT(numLats))
798	iret = NF90_GET_VAR(fid, latid, LAT(:), start = (/ 1, 1, 1 /), count = (/ 1, numLats, 1 /))
799
800	user_dlon = (LON(numLons) - LON(1)) / (numLons - 1)
801	user_dlat = (LAT(numLats) - LAT(1)) / (numLats - 1)
802
803	DEALLOCATE(LON,LAT)
804
805	ENDIF
806	! Unable to know from where to get the information
807	user_known_x = 1
808	user_known_y = 1
809	!
810	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "corner_lats", corner_lats)
811	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "corner_lons", corner_lons)
812	user_known_lat = corner_lats(1)
813	user_known_lon = corner_lons(1)
814	!
815	! Read mapfactor, land mask and orography
816	!
817	!
818	! Allocation
819	!
820	ALLOCATE(mapfac_x(iim,jjm))
821	ALLOCATE(mapfac_y(iim,jjm))
822	ALLOCATE(dxwrf(iim,jjm))
823	ALLOCATE(dywrf(iim,jjm))
824	!
825	varname = "MAPFAC_MX"
826	iret = NF90_INQ_VARID (fid, varname, iv)
827	IF (iret /= NF90_NOERR) THEN
828	CALL ipslerr (3,'WRFdomain_Read',"Could not find variable ", varname," ")
829	ELSE
830	iret = NF90_GET_VAR (fid,iv,mapfac_x)
831	ENDIF
832	varname = "MAPFAC_MY"
833	iret = NF90_INQ_VARID (fid, varname, iv)
834	IF (iret /= NF90_NOERR) THEN
835	CALL ipslerr (3,'WRFdomain_Read',"Could not find variable ", varname," ")
836	ELSE
837	iret = NF90_GET_VAR (fid,iv,mapfac_y)
838	ENDIF
839	!
840	! Initilize the projection
841	!
842	if (current_proj == PROJ_LATLON) then
843	call map_set(current_proj, proj_stack(idom), &
844	lat1=user_known_lat, &
845	lon1=user_known_lon, &
846	knowni=user_known_x, &
847	knownj=user_known_y, &
848	latinc=user_dlat, &
849	loninc=user_dlon, &
850	r_earth=R_Earth)
851
852	else if (current_proj == PROJ_MERC) then
853	call map_set(current_proj, proj_stack(idom), &
854	truelat1=user_truelat1, &
855	lat1=user_known_lat, &
856	lon1=user_known_lon, &
857	knowni=user_known_x, &
858	knownj=user_known_y, &
859	dx=user_dxkm, &
860	r_earth=R_Earth)
861
862	else if (current_proj == PROJ_CYL) then
863	call ipslerr(3,"grid_initproj",'Should not have PROJ_CYL as projection for',&
864	'source data in push_source_projection()', " ")
865
866	else if (current_proj == PROJ_CASSINI) then
867	call ipslerr(3,"grid_initproj",'Should not have PROJ_CASSINI as projection for', &
868	'source data in push_source_projection()', " ")
869
870	else if (current_proj == PROJ_LC) then
871	call map_set(current_proj, proj_stack(idom), &
872	truelat1=user_truelat1, &
873	truelat2=user_truelat2, &
874	stdlon=user_stand_lon, &
875	lat1=user_known_lat, &
876	lon1=user_known_lon, &
877	knowni=user_known_x, &
878	knownj=user_known_y, &
879	dx=user_dxkm, &
880	r_earth=R_Earth)
881
882	else if (current_proj == PROJ_ALBERS_NAD83) then
883	call map_set(current_proj, proj_stack(idom), &
884	truelat1=user_truelat1, &
885	truelat2=user_truelat2, &
886	stdlon=user_stand_lon, &
887	lat1=user_known_lat, &
888	lon1=user_known_lon, &
889	knowni=user_known_x, &
890	knownj=user_known_y, &
891	dx=user_dxkm, &
892	r_earth=R_Earth)
893
894	else if (current_proj == PROJ_PS) then
895	call map_set(current_proj, proj_stack(idom), &
896	truelat1=user_truelat1, &
897	stdlon=user_stand_lon, &
898	lat1=user_known_lat, &
899	lon1=user_known_lon, &
900	knowni=user_known_x, &
901	knownj=user_known_y, &
902	dx=user_dxkm, &
903	r_earth=R_Earth)
904
905	else if (current_proj == PROJ_PS_WGS84) then
906	call map_set(current_proj, proj_stack(idom), &
907	truelat1=user_truelat1, &
908	stdlon=user_stand_lon, &
909	lat1=user_known_lat, &
910	lon1=user_known_lon, &
911	knowni=user_known_x, &
912	knownj=user_known_y, &
913	dx=user_dxkm, &
914	r_earth=R_Earth)
915
916	else if (current_proj == PROJ_GAUSS) then
917	call map_set(current_proj, proj_stack(idom), &
918	lat1=user_known_lat, &
919	lon1=user_known_lon, &
920	nlat=nint(user_dlat), &
921	loninc=user_dlon, &
922	r_earth=R_Earth)
923
924	else if (current_proj == PROJ_ROTLL) then
925	call ipslerr(3 ,"grid_initproj",'Should not have PROJ_ROTLL as projection for', &
926	'source data in push_source_projection() as not yet implemented', '')
927	end if
928	!
929	! Transform the mapfactors into dx and dy to be used for the description of the polygons and
930	! interpolations.
931	!
932	DO i=1,iim
933	DO j=1,jjm
934	!
935	IF (proj_stack(idom)%code /= PROJ_LATLON ) THEN
936	dx = proj_stack(idom)%dx
937	! Some projections in WRF do not store dy, in that case dy=dx.
938	IF ( proj_stack(idom)%dy > 0 ) THEN
939	dy = proj_stack(idom)%dy
940	ELSE
941	dy = proj_stack(idom)%dx
942	ENDIF
943	dxwrf(i,j) = dx/mapfac_x(i,j)
944	dywrf(i,j) = dy/mapfac_y(i,j)
945	ELSE
946	!
947	! The LatLon projection is also a special case as here it is not the dx and dy
948	! which are stored in the projection file but the increments in Lon and Lat.
949	!
950	dtx = proj_stack(idom)%loninc
951	dty = proj_stack(idom)%latinc
952	coslat = COS(lat(j) * pi/180. )
953	dxwrf(i,j) = dtx * pi/180. * R_Earth * coslat
954	dywrf(i,j) = dty * pi/180. * R_Earth
955	!
956	ENDIF
957	!
958	ENDDO
959	ENDDO
960	!
961	END SUBROUTINE grid_initproj
962	!
963	!
964	!
965	!=========================================================================================
966	!
967	SUBROUTINE grid_tolola_scal (ri, rj, lon, lat)
968	!
969	!
970	! Argument
971	REAL(r_std), INTENT(in) :: ri, rj
972	REAL(r_std), INTENT(out) :: lon, lat
973	!
974	!
975	IF ( proj_stack(1)%code < undef_int ) THEN
976	!
977	CALL ij_to_latlon(proj_stack(1), ri, rj, lat, lon)
978	!
979	ELSE
980	CALL ipslerr(3, "grid_tolola_scal", "Projection not initilized"," "," ")
981	ENDIF
982	!
983	END SUBROUTINE grid_tolola_scal
984	!
985	!=========================================================================================
986	!
987	SUBROUTINE grid_tolola_1d (ri, rj, lon, lat)
988	!
989	!
990	! Argument
991	REAL(r_std), INTENT(in), DIMENSION(:) :: ri, rj
992	REAL(r_std), INTENT(out), DIMENSION(:) :: lon, lat
993	!
994	! Local
995	INTEGER :: i, imax
996	!
997	imax=SIZE(lon)
998	!
999	IF ( proj_stack(1)%code < undef_int ) THEN
1000	DO i=1,imax
1001	!
1002	CALL ij_to_latlon(proj_stack(1), ri(i), rj(i), lat(i), lon(i))
1003	!
1004	ENDDO
1005	ELSE
1006	CALL ipslerr(3, "grid_tolola_1d", "Projection not initilized"," "," ")
1007	ENDIF
1008	!
1009	END SUBROUTINE grid_tolola_1d
1010	!
1011	!=========================================================================================
1012	!
1013	SUBROUTINE grid_tolola_2d (ri, rj, lon, lat)
1014	!
1015	!
1016	! Argument
1017	REAL(r_std), INTENT(in), DIMENSION(:,:) :: ri, rj
1018	REAL(r_std), INTENT(out), DIMENSION(:,:) :: lon, lat
1019	!
1020	! Local
1021	INTEGER :: i, imax, j, jmax
1022	!
1023	imax=SIZE(lon,DIM=1)
1024	jmax=SIZE(lon,DIM=2)
1025	!
1026	IF ( proj_stack(1)%code < undef_int ) THEN
1027	DO i=1,imax
1028	DO j=1,jmax
1029	!
1030	CALL ij_to_latlon(proj_stack(1), ri(i,j), rj(i,j), lat(i,j), lon(i,j))
1031	!
1032	ENDDO
1033	ENDDO
1034	ELSE
1035	CALL ipslerr(3, "grid_tolola_2d", "Projection not initilized"," "," ")
1036	ENDIF
1037	!
1038	END SUBROUTINE grid_tolola_2d
1039	!
1040	!=========================================================================================
1041	!
1042	SUBROUTINE grid_toij_scal (lon, lat, ri, rj)
1043	!
1044	!
1045	! Argument
1046	REAL(r_std), INTENT(in) :: lon, lat
1047	REAL(r_std), INTENT(out) :: ri, rj
1048	!
1049	!
1050	IF ( proj_stack(1)%code < undef_int ) THEN
1051	!
1052	CALL latlon_to_ij(proj_stack(1), lat, lon, ri, rj)
1053	!
1054	ELSE
1055	CALL ipslerr(3, "grid_toij_scal", "Projection not initilized"," "," ")
1056	ENDIF
1057	!
1058	END SUBROUTINE grid_toij_scal
1059	!
1060	!=========================================================================================
1061	!
1062	SUBROUTINE grid_toij_1d (lon, lat, ri, rj)
1063	!
1064	!
1065	! Argument
1066	REAL(r_std), INTENT(in), DIMENSION(:) :: lon, lat
1067	REAL(r_std), INTENT(out), DIMENSION(:) :: ri, rj
1068	!
1069	! Local
1070	INTEGER :: i, imax
1071	!
1072	imax=SIZE(lon)
1073	!
1074	IF ( proj_stack(1)%code < undef_int ) THEN
1075	DO i=1,imax
1076	!
1077	CALL latlon_to_ij(proj_stack(1), lat(i), lon(i), ri(i), rj(i))
1078	!
1079	ENDDO
1080	ELSE
1081	CALL ipslerr(3, "grid_toij_1d", "Projection not initilized"," "," ")
1082	ENDIF
1083	!
1084	END SUBROUTINE grid_toij_1d
1085	!
1086	!=========================================================================================
1087	!
1088	SUBROUTINE grid_toij_2d (lon, lat, ri, rj)
1089	!
1090	!
1091	! Argument
1092	REAL(r_std), INTENT(in), DIMENSION(:,:) :: lon, lat
1093	REAL(r_std), INTENT(out), DIMENSION(:,:) :: ri, rj
1094	!
1095	! Local
1096	INTEGER :: i, imax, j, jmax
1097	!
1098	imax=SIZE(lon,DIM=1)
1099	jmax=SIZE(lon,DIM=2)
1100	!
1101	IF ( proj_stack(1)%code < undef_int ) THEN
1102	DO i=1,imax
1103	DO j=1,jmax
1104	!
1105	CALL latlon_to_ij(proj_stack(1), lat(i,j), lon(i,j), ri(i,j), rj(i,j))
1106	!
1107	ENDDO
1108	ENDDO
1109	ELSE
1110	CALL ipslerr(3, "grid_toij_2d", "Projection not initilized"," "," ")
1111	ENDIF
1112	!
1113	END SUBROUTINE grid_toij_2d
1114	!
1115	!
1116	!=========================================================================================
1117	!
1118	!
1119	END MODULE grid

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

Download in other formats: