Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_global/grid.f90 @ 6393

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

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

Download in other formats: