Context Navigation

source: branches/ORCHIDEE_3_CMIP6/ORCHIDEE/src_global/grid.f90 @ 8765

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

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

Download in other formats: