Context Navigation

source: tags/ORCHIDEE_2_1/ORCHIDEE/src_global/grid.f90

Last change on this file was 5559, checked in by josefine.ghattas, 6 years ago
Enhancement on grid module: merge variables grid_type and GridType? meaning the same. Only grid_type is used now, this one is choosen because it is the same as used in LMDZ.
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 37.7 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	ALLOCATE(neighbours_g(nbp_glo,NbNeighb))
364	ALLOCATE(headings_g(nbp_glo,NbNeighb))
365	ALLOCATE(seglength_g(nbp_glo,NbSegments))
366	ALLOCATE(corners_g(nbp_glo,NbSegments,2))
367	ALLOCATE(area_g(nbp_glo))
368	!
369	! TEMPORARY
370	!
371	ALLOCATE(resolution_g(nbp_glo,2))
372	!
373	! Allocate other variables
374	!
375	ALLOCATE(lalo_g(nbp_glo,2), contfrac_g(nbp_glo),index_g(nbp_glo))
376	ALLOCATE(lon_g(iim_g, jjm_g), lat_g(iim_g, jjm_g))
377	!
378	END SUBROUTINE grid_allocate_glo
379	!!
380	!! =============================================================================================================================
381	!! SUBROUTINE: grid_stuff
382	!!
383	!>\BRIEF transfers the global horizontal grid information to ORCHIDEE in the case of grid regular in Longitude
384	!! and Latitude.
385	!!
386	!! DESCRIPTION:
387	!!
388	!!
389	!! This subroutine is called by intersurf_main_2d or any driver of the model.
390	!!
391	!! \n
392	!_ ==============================================================================================================================
393	!!
394	SUBROUTINE grid_stuff (npts_glo, iim, jjm, grid_lon, grid_lat, kindex, contfrac_tmp)
395	!
396	! 0 interface
397	!
398	IMPLICIT NONE
399	!
400	! 0.1 input !
401
402	! Domain size
403	INTEGER(i_std), INTENT(in) :: npts_glo
404	! Size of cartesian grid
405	INTEGER(i_std), INTENT(in) :: iim, jjm
406	! Longitudes on cartesian grid
407	REAL(r_std), DIMENSION(iim,jjm), INTENT(in) :: grid_lon
408	! Latitudes on cartesian grid
409	REAL(r_std), DIMENSION(iim,jjm), INTENT(in) :: grid_lat
410	! Index of land point on 2D map (in local position)
411	INTEGER(i_std), DIMENSION(npts_glo), INTENT(in) :: kindex
412	! The fraction of continent in the grid box [0-1]
413	REAL(r_std), DIMENSION(npts_glo), OPTIONAL, INTENT(in) :: contfrac_tmp
414	!
415	!
416	! =========================================================================
417
418	IF ( printlev >= 4 ) WRITE(numout,*) 'Entering grid_stuff'
419
420	! default resolution
421	IF ( printlev >=2 ) WRITE(numout,*) 'grid stuff: default resolution (m): ',default_resolution
422	!
423	!-
424	IF (is_root_prc) THEN
425	!
426	CALL grid_topolylist(NbSegments, npts_glo, iim, jjm, grid_lon, grid_lat, kindex, &
427	& global, corners_g, neighbours_g, headings_g, seglength_g, area_g, ilandindex, jlandindex)
428	!
429	IF (PRESENT(contfrac_tmp)) THEN
430	!
431	! Transfer the contfrac into the array managed in this module.
432	!
433	contfrac_g(:) = contfrac_tmp(:)
434	ENDIF
435	!
436	ENDIF
437	!
438	! With this the description of the grid is complete and the information
439	! can be scattered to all processors.
440	!
441	CALL grid_scatter()
442	!
443	IF ( printlev >= 3 ) WRITE(numout,*) 'Leaving grid_stuff'
444
445	END SUBROUTINE grid_stuff
446	!!
447	!! =============================================================================================================================
448	!! SUBROUTINE: grid_topolylist
449	!!
450	!>\BRIEF This routine transforms a regular grid into a list of polygons which are defined by the following
451	!! quantities :
452	!!
453	!! corners : the n vertices of the polugon in longitude and latitude
454	!! neighbours : the neighbouring land grid box for each of the vertices and segments
455	!! headings : the direction in which the neighbour is
456	!! seglength : the lenght of each segment
457	!! area : the area of the polygon
458	!! ilindex, jlindex : provides the i,j coordinates of the mesh in the global grid.
459	!!
460	!! DESCRIPTION:
461	!!
462	!! \n
463	!_ ==============================================================================================================================
464	!!
465	SUBROUTINE grid_topolylist(nbseg, nland, iim, jjm, grid_lon, grid_lat, kindex, &
466	& globalg, corners_loc, neighbours_loc, headings_loc, seglength_loc, &
467	& area_loc, ilindex_loc, jlindex_loc)
468	!
469	! 0 interface
470	!
471	IMPLICIT NONE
472	!
473	! 0.1 input !
474	! Number of segments for each polygon
475	INTEGER(i_std), INTENT(in) :: nbseg
476	! Number of land points on the grid
477	INTEGER(i_std), INTENT(in) :: nland
478	! Size of cartesian grid
479	INTEGER(i_std), INTENT(in) :: iim, jjm
480	! Longitudes on cartesian grid
481	REAL(r_std), DIMENSION(iim,jjm), INTENT(in) :: grid_lon
482	! Latitudes on cartesian grid
483	REAL(r_std), DIMENSION(iim,jjm), INTENT(in) :: grid_lat
484	! Index of land point on 2D map (in local position)
485	INTEGER(i_std), DIMENSION(nland), INTENT(in) :: kindex
486	!
487	! 0.2 Output
488	!
489	LOGICAL, INTENT(inout) :: globalg
490	!
491	REAL(r_std), DIMENSION(nland,nbseg,2), INTENT(out) :: corners_loc
492	INTEGER(i_std), DIMENSION(nland,nbseg*2), INTENT(out) :: neighbours_loc
493	REAL(r_std), DIMENSION(nland,nbseg*2), INTENT(out) :: headings_loc
494	REAL(r_std), DIMENSION(nland,nbseg), INTENT(out) :: seglength_loc
495	REAL(r_std), DIMENSION(nland), INTENT(out) :: area_loc
496	INTEGER(i_std), DIMENSION(nland), INTENT(out) :: ilindex_loc, jlindex_loc
497	!
498	! 0.3 Local variables
499	!
500	INTEGER(i_std) :: i, is, iss
501	REAL(r_std), DIMENSION(nland,2) :: center
502	REAL(r_std) :: maxdellon, mindellon, maxlon, minlon
503	REAL(r_std), DIMENSION(nland,nbseg*2) :: lonpoly, latpoly
504	!
505	IF ( grid_type == regular_lonlat ) THEN
506	!
507	! If we are in regular Lon Lat, then we test just the longitude and see if we span 0-360deg.
508	!
509	maxdellon=MAXVAL(ABS(grid_lon(1:iim-1,1)-grid_lon(2:iim,1)))
510	mindellon=MINVAL(ABS(grid_lon(1:iim-1,1)-grid_lon(2:iim,1)))
511	maxlon=MAXVAL(grid_lon(1:iim,1))
512	minlon=MINVAL(grid_lon(1:iim,1))
513	!
514	! test if it could be a global grid on 0 -> 360
515	!
516	IF ( minlon > 0 .AND. maxlon > 180 ) THEN
517	IF ( (minlon - (maxdellon + min_sechiba)) <= 0 .AND. (maxlon + (maxdellon + min_sechiba)) >= 360) THEN
518	globalg = .TRUE.
519	ELSE
520	globalg = .FALSE.
521	ENDIF
522	!
523	! Test if it could be a -180 to 180 grid
524	!
525	ELSE IF ( minlon < 0 .AND. maxlon > 0 ) THEN
526	IF ( (minlon - (maxdellon + min_sechiba)) <= -180 .AND. (maxlon + (maxdellon + min_sechiba)) >= 180) THEN
527	globalg = .TRUE.
528	ELSE
529	globalg = .FALSE.
530	ENDIF
531	!
532	! If neither condition is met then it cannot be global.
533	!
534	ELSE
535	globalg = .FALSE.
536	ENDIF
537	ELSE IF ( grid_type == regular_xy ) THEN
538	!
539	! The hypothesis is that if we are using grid regular_xy then it is not a global grid
540	!
541	globalg = .FALSE.
542	ELSE
543	STOP "Unknown grid"
544	ENDIF
545	!
546	! 2.0 Transform the grid into a list of polygones while keeping the neighbour relations
547	! between these polygones.
548	!
549	! Each polygone starts with a vertex and alternates vertices and mid-points of segments.
550	!
551	IF (nland == 1) THEN
552	CALL haversine_singlepointploy(iim, jjm, grid_lon, grid_lat, nland, kindex, globalg, &
553	& nbseg, lonpoly, latpoly, center, &
554	& neighbours_loc, ilindex_loc, jlindex_loc)
555	ELSE IF ( grid_type == regular_lonlat ) THEN
556	CALL haversine_reglatlontoploy(iim, jjm, grid_lon, grid_lat, nland, kindex, globalg, &
557	& nbseg, lonpoly, latpoly, center, &
558	& neighbours_loc, ilindex_loc, jlindex_loc)
559	ELSE IF ( grid_type == regular_xy ) THEN
560	CALL haversine_regxytoploy(iim, jjm, grid_lon, grid_lat, nland, kindex, proj_stack, &
561	& nbseg, lonpoly, latpoly, center, &
562	& neighbours_loc, ilindex_loc, jlindex_loc)
563	ELSE
564	STOP "Unknown grid"
565	ENDIF
566	!
567	! Save the longitude and latitudes nbseg corners (=vertices) of the polygones
568	!
569	DO i=1,nland
570	DO is=1,nbseg
571	iss=(is-1)*2+1
572	corners_loc(i,is,1) = lonpoly(i,iss)
573	corners_loc(i,is,2) = latpoly(i,iss)
574	ENDDO
575	ENDDO
576	!
577	! Get the heading normal to the 4 segments and through the 4 corners.
578	!
579	CALL haversine_polyheadings(nland, nbseg, lonpoly, latpoly, center, headings_loc)
580	!
581	! Order the points of the polygone in clockwise order Starting with the northern most
582	!
583	CALL haversine_polysort(nland, nbseg, lonpoly, latpoly, headings_loc, neighbours_loc)
584	!
585	! Compute the segment length and area.
586	! For the RegLonLat we have specific calculations for seglength and area.
587	! For projected regular grids we use the great cicle assumption for the segments
588	! but the projected area.
589	! For unstructured grid we use the most general routines.
590	!
591	IF ( grid_type == regular_lonlat ) THEN
592	CALL haversine_laloseglen(nland, nbseg, lonpoly, latpoly, seglength_loc)
593	CALL haversine_laloarea(nland, nbseg, seglength_loc, area_loc)
594	ELSE IF ( grid_type == regular_xy ) THEN
595	CALL haversine_polyseglen(nland, nbseg, lonpoly, latpoly, seglength_loc)
596	CALL haversine_xyarea(nland, nbseg, ilindex_loc, jlindex_loc, dxwrf, dywrf, area_loc)
597	ELSE
598	CALL haversine_polyseglen(nland, nbseg, lonpoly, latpoly, seglength_loc)
599	CALL haversine_polyarea(nland, nbseg, lonpoly, latpoly, area_loc)
600	ENDIF
601
602	END SUBROUTINE grid_topolylist
603	!!
604	!!
605	!!
606	!! =============================================================================================================================
607	!! SUBROUTINE: grid_scatter
608	!!
609	!>\BRIEF Scatter the grid information so that each processor knows the characteristics of the grid it works on.
610	!!
611	!! DESCRIPTION:
612	!!
613	!!
614	!! The grid information has been computed for the entire grid on the root processor. Now we give each processor
615	!! the information of the piece of the grid it works on. This concerns the following variables describing the grid :
616	!! - area
617	!! - resolution
618	!! - neighbours
619	!! - contfrac : fraction of continent
620	!!
621	!! Should ilandindex and jlandindex not b initialized, we catch-up here. This field is the same on all processors.
622	!!
623	!! TODO :
624	!! This code should get the grid describing fields as arguments and then writem into the *_g variables on
625	!! root_prc before scattering. This would allow to compute the grid characteristics in any subroutine
626	!! fore calling grid_scatter.
627	!!
628	!!
629	!!
630	!! \n
631	!_ ==============================================================================================================================
632	!!
633	!!
634	SUBROUTINE grid_scatter()
635	!
636	!
637	INTEGER(i_std) :: i, ip, jp
638	!
639	IF ( MAXVAL(ilandindex) < 0 .AND. MAXVAL(jlandindex) < 0 ) THEN
640	DO i = 1, nbp_glo
641	!
642	! 1 find numbers of the latitude and longitude of each point
643	!
644
645	! index of latitude
646	jp = INT( (index_g(i)-1) /iim_g ) + 1
647
648	! index of longitude
649	ip = index_g(i) - ( jp-1 ) * iim_g
650	!
651	! Save this information for usage in other modules.
652	!
653	ilandindex(i)=ip
654	jlandindex(i)=jp
655	!
656	ENDDO
657	ENDIF
658	!
659	CALL scatter(neighbours_g, neighbours)
660	CALL scatter(contfrac_g, contfrac)
661	CALL scatter(headings_g, headings)
662	CALL scatter(seglength_g, seglength)
663	CALL scatter(corners_g, corners)
664	CALL scatter(area_g, area)
665	!
666	! TEMPORARY section for resolution
667	!
668	IF ( is_root_prc) THEN
669	IF ( grid_type == regular_lonlat .OR. grid_type == regular_xy ) THEN
670	resolution_g(:,1) = (seglength_g(:,1)+seglength_g(:,3))/2.0
671	resolution_g(:,2) = (seglength_g(:,2)+seglength_g(:,4))/2.0
672	ELSE
673	CALL ipslerr(3, "grid_scatter", "unsupported grid type.",&
674	& "grid_scatter can only be called for regular grids.",&
675	& "")
676	ENDIF
677	ENDIF
678	CALL scatter(resolution_g, resolution)
679
680	! Copy variable global from root processor to all prossesors
681	CALL bcast(global)
682
683	!
684	!
685	IF ( printlev >=4 ) THEN
686	WRITE(numout,*) 'grid_scatter > seglength = ', seglength(1,:)
687	WRITE(numout,*) 'grid_scatter > neighbours = ', neighbours(1,:)
688	WRITE(numout,*) 'grid_scatter > contfrac = ', contfrac(1)
689	WRITE(numout,*) 'grid_scatter > area = ', area(1)
690	WRITE(numout,*) 'grid_scatter > global = ', global
691	ENDIF
692	!
693	END SUBROUTINE grid_scatter
694	!!
695	!!
696	!! =============================================================================================================================
697	!! SUBROUTINE: grid_initproj
698	!!
699	!>\BRIEF Routine to initialise the projection
700	!!
701	!! DESCRIPTION:
702	!!
703	!!
704	!! This subroutine is called by the routine whichs ets-up th grid on which ORCHIDEE is to run.
705	!! The aim is to set-upu the projection so that all the grid variables needed by ORCHIDEE can
706	!! be computed in grid_stuff_regxy
707	!!
708	!! \n
709	!_ ==============================================================================================================================
710	!!
711	!!
712	SUBROUTINE grid_initproj (fid, iim, jjm)
713	!
714	!
715	! 0 interface
716	!
717	IMPLICIT NONE
718	!
719	! 0.1 input !
720	!
721	! Domain size
722	INTEGER(i_std), INTENT(in) :: fid
723	INTEGER(i_std), INTENT(in) :: iim, jjm
724	!
725	! 0.2 Local variables
726	!
727	INTEGER(i_std) :: current_proj, idom, iret, lonid, latid, numLons, numLats
728	INTEGER, DIMENSION(nf90_max_var_dims) :: dimIDs
729	REAL(r_std) :: user_stand_lon, user_truelat1, user_truelat2, user_dxkm, user_dykm
730	REAL(r_std) :: user_dlat, user_dlon, user_known_x, user_known_y, user_known_lat, user_known_lon
731	REAL(r_std), DIMENSION(16) :: corner_lons, corner_lats
732	!
733	INTEGER(i_std) :: iv, i, j
734	CHARACTER(LEN=20) :: varname
735	REAL(r_std) :: dx, dy, dtx, dty, coslat
736	REAL(r_std), ALLOCATABLE, DIMENSION (:) :: LON, LAT
737	REAL(r_std), ALLOCATABLE, DIMENSION (:,:) :: mapfac_x, mapfac_y
738	!
739	!
740	! Only one domain is possible for the moment
741	!
742	idom=1
743	CALL map_init(proj_stack(idom))
744	!
745	! Does ORCHIDEE have the same Earth Radius as the map projection ?
746	!
747	IF ( ABS(R_Earth-EARTH_RADIUS_M) > 0.1 ) THEN
748	WRITE(,) "Earth Radius in WRF : ", EARTH_RADIUS_M
749	WRITE(,) "Earth Radius in ORCHIDEE : ", R_Earth
750	CALL ipslerr (3,'grid_initproj','The Earth radius is not the same in the projection module and ORCHIDEE',&
751	& " ", " ")
752	ENDIF
753	!
754	! Get parameters of the projection from the netCDF file
755	!
756	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "MAP_PROJ", current_proj)
757	!
758	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "STAND_LON", user_stand_lon)
759	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "TRUELAT1", user_truelat1)
760	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "TRUELAT2", user_truelat2)
761	!
762	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "DX", user_dxkm)
763	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "DY", user_dykm)
764	user_dlat = undef
765	user_dlon = undef
766	!
767	IF ( current_proj == PROJ_LATLON ) THEN
768	!
769	iret = NF90_inq_VARID(fid, "XLONG_M",lonid)
770	iret = NF90_INQUIRE_VARIABLE(fid, lonid, dimids = dimIDs)
771	iret = NF90_INQUIRE_DIMENSION(fid, dimIDs(1), len = numLons)
772	iret = NF90_INQUIRE_DIMENSION(fid, dimIDs(2), len = numLats)
773	ALLOCATE(LON(numLons))
774	iret = NF90_GET_VAR(fid, lonid, LON(:), start = (/ 1, 1, 1 /), count = (/ numLons, 1, 1 /))
775
776	iret = NF90_inq_VARID(fid, "XLAT_M",latid)
777	ALLOCATE(LAT(numLats))
778	iret = NF90_GET_VAR(fid, latid, LAT(:), start = (/ 1, 1, 1 /), count = (/ 1, numLats, 1 /))
779
780	user_dlon = (LON(numLons) - LON(1)) / (numLons - 1)
781	user_dlat = (LAT(numLats) - LAT(1)) / (numLats - 1)
782
783	DEALLOCATE(LON,LAT)
784
785	ENDIF
786	! Unable to know from where to get the information
787	user_known_x = 1
788	user_known_y = 1
789	!
790	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "corner_lats", corner_lats)
791	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "corner_lons", corner_lons)
792	user_known_lat = corner_lats(1)
793	user_known_lon = corner_lons(1)
794	!
795	! Read mapfactor, land mask and orography
796	!
797	!
798	! Allocation
799	!
800	ALLOCATE(mapfac_x(iim,jjm))
801	ALLOCATE(mapfac_y(iim,jjm))
802	ALLOCATE(dxwrf(iim,jjm))
803	ALLOCATE(dywrf(iim,jjm))
804	!
805	varname = "MAPFAC_MX"
806	iret = NF90_INQ_VARID (fid, varname, iv)
807	IF (iret /= NF90_NOERR) THEN
808	CALL ipslerr (3,'WRFdomain_Read',"Could not find variable ", varname," ")
809	ELSE
810	iret = NF90_GET_VAR (fid,iv,mapfac_x)
811	ENDIF
812	varname = "MAPFAC_MY"
813	iret = NF90_INQ_VARID (fid, varname, iv)
814	IF (iret /= NF90_NOERR) THEN
815	CALL ipslerr (3,'WRFdomain_Read',"Could not find variable ", varname," ")
816	ELSE
817	iret = NF90_GET_VAR (fid,iv,mapfac_y)
818	ENDIF
819	!
820	! Initilize the projection
821	!
822	if (current_proj == PROJ_LATLON) then
823	call map_set(current_proj, proj_stack(idom), &
824	lat1=user_known_lat, &
825	lon1=user_known_lon, &
826	knowni=user_known_x, &
827	knownj=user_known_y, &
828	latinc=user_dlat, &
829	loninc=user_dlon, &
830	r_earth=R_Earth)
831
832	else if (current_proj == PROJ_MERC) then
833	call map_set(current_proj, proj_stack(idom), &
834	truelat1=user_truelat1, &
835	lat1=user_known_lat, &
836	lon1=user_known_lon, &
837	knowni=user_known_x, &
838	knownj=user_known_y, &
839	dx=user_dxkm, &
840	r_earth=R_Earth)
841
842	else if (current_proj == PROJ_CYL) then
843	call ipslerr(3,"grid_initproj",'Should not have PROJ_CYL as projection for',&
844	'source data in push_source_projection()', " ")
845
846	else if (current_proj == PROJ_CASSINI) then
847	call ipslerr(3,"grid_initproj",'Should not have PROJ_CASSINI as projection for', &
848	'source data in push_source_projection()', " ")
849
850	else if (current_proj == PROJ_LC) then
851	call map_set(current_proj, proj_stack(idom), &
852	truelat1=user_truelat1, &
853	truelat2=user_truelat2, &
854	stdlon=user_stand_lon, &
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_ALBERS_NAD83) then
863	call map_set(current_proj, proj_stack(idom), &
864	truelat1=user_truelat1, &
865	truelat2=user_truelat2, &
866	stdlon=user_stand_lon, &
867	lat1=user_known_lat, &
868	lon1=user_known_lon, &
869	knowni=user_known_x, &
870	knownj=user_known_y, &
871	dx=user_dxkm, &
872	r_earth=R_Earth)
873
874	else if (current_proj == PROJ_PS) then
875	call map_set(current_proj, proj_stack(idom), &
876	truelat1=user_truelat1, &
877	stdlon=user_stand_lon, &
878	lat1=user_known_lat, &
879	lon1=user_known_lon, &
880	knowni=user_known_x, &
881	knownj=user_known_y, &
882	dx=user_dxkm, &
883	r_earth=R_Earth)
884
885	else if (current_proj == PROJ_PS_WGS84) then
886	call map_set(current_proj, proj_stack(idom), &
887	truelat1=user_truelat1, &
888	stdlon=user_stand_lon, &
889	lat1=user_known_lat, &
890	lon1=user_known_lon, &
891	knowni=user_known_x, &
892	knownj=user_known_y, &
893	dx=user_dxkm, &
894	r_earth=R_Earth)
895
896	else if (current_proj == PROJ_GAUSS) then
897	call map_set(current_proj, proj_stack(idom), &
898	lat1=user_known_lat, &
899	lon1=user_known_lon, &
900	nlat=nint(user_dlat), &
901	loninc=user_dlon, &
902	r_earth=R_Earth)
903
904	else if (current_proj == PROJ_ROTLL) then
905	call ipslerr(3 ,"grid_initproj",'Should not have PROJ_ROTLL as projection for', &
906	'source data in push_source_projection() as not yet implemented', '')
907	end if
908	!
909	! Transform the mapfactors into dx and dy to be used for the description of the polygons and
910	! interpolations.
911	!
912	DO i=1,iim
913	DO j=1,jjm
914	!
915	IF (proj_stack(idom)%code /= PROJ_LATLON ) THEN
916	dx = proj_stack(idom)%dx
917	! Some projections in WRF do not store dy, in that case dy=dx.
918	IF ( proj_stack(idom)%dy > 0 ) THEN
919	dy = proj_stack(idom)%dy
920	ELSE
921	dy = proj_stack(idom)%dx
922	ENDIF
923	dxwrf(i,j) = dx/mapfac_x(i,j)
924	dywrf(i,j) = dy/mapfac_y(i,j)
925	ELSE
926	!
927	! The LatLon projection is also a special case as here it is not the dx and dy
928	! which are stored in the projection file but the increments in Lon and Lat.
929	!
930	dtx = proj_stack(idom)%loninc
931	dty = proj_stack(idom)%latinc
932	coslat = COS(lat(j) * pi/180. )
933	dxwrf(i,j) = dtx * pi/180. * R_Earth * coslat
934	dywrf(i,j) = dty * pi/180. * R_Earth
935	!
936	ENDIF
937	!
938	ENDDO
939	ENDDO
940	!
941	END SUBROUTINE grid_initproj
942	!
943	!
944	!
945	!=========================================================================================
946	!
947	SUBROUTINE grid_tolola_scal (ri, rj, lon, lat)
948	!
949	!
950	! Argument
951	REAL(r_std), INTENT(in) :: ri, rj
952	REAL(r_std), INTENT(out) :: lon, lat
953	!
954	!
955	IF ( proj_stack(1)%code < undef_int ) THEN
956	!
957	CALL ij_to_latlon(proj_stack(1), ri, rj, lat, lon)
958	!
959	ELSE
960	CALL ipslerr(3, "grid_tolola_scal", "Projection not initilized"," "," ")
961	ENDIF
962	!
963	END SUBROUTINE grid_tolola_scal
964	!
965	!=========================================================================================
966	!
967	SUBROUTINE grid_tolola_1d (ri, rj, lon, lat)
968	!
969	!
970	! Argument
971	REAL(r_std), INTENT(in), DIMENSION(:) :: ri, rj
972	REAL(r_std), INTENT(out), DIMENSION(:) :: lon, lat
973	!
974	! Local
975	INTEGER :: i, imax
976	!
977	imax=SIZE(lon)
978	!
979	IF ( proj_stack(1)%code < undef_int ) THEN
980	DO i=1,imax
981	!
982	CALL ij_to_latlon(proj_stack(1), ri(i), rj(i), lat(i), lon(i))
983	!
984	ENDDO
985	ELSE
986	CALL ipslerr(3, "grid_tolola_1d", "Projection not initilized"," "," ")
987	ENDIF
988	!
989	END SUBROUTINE grid_tolola_1d
990	!
991	!=========================================================================================
992	!
993	SUBROUTINE grid_tolola_2d (ri, rj, lon, lat)
994	!
995	!
996	! Argument
997	REAL(r_std), INTENT(in), DIMENSION(:,:) :: ri, rj
998	REAL(r_std), INTENT(out), DIMENSION(:,:) :: lon, lat
999	!
1000	! Local
1001	INTEGER :: i, imax, j, jmax
1002	!
1003	imax=SIZE(lon,DIM=1)
1004	jmax=SIZE(lon,DIM=2)
1005	!
1006	IF ( proj_stack(1)%code < undef_int ) THEN
1007	DO i=1,imax
1008	DO j=1,jmax
1009	!
1010	CALL ij_to_latlon(proj_stack(1), ri(i,j), rj(i,j), lat(i,j), lon(i,j))
1011	!
1012	ENDDO
1013	ENDDO
1014	ELSE
1015	CALL ipslerr(3, "grid_tolola_2d", "Projection not initilized"," "," ")
1016	ENDIF
1017	!
1018	END SUBROUTINE grid_tolola_2d
1019	!
1020	!=========================================================================================
1021	!
1022	SUBROUTINE grid_toij_scal (lon, lat, ri, rj)
1023	!
1024	!
1025	! Argument
1026	REAL(r_std), INTENT(in) :: lon, lat
1027	REAL(r_std), INTENT(out) :: ri, rj
1028	!
1029	!
1030	IF ( proj_stack(1)%code < undef_int ) THEN
1031	!
1032	CALL latlon_to_ij(proj_stack(1), lat, lon, ri, rj)
1033	!
1034	ELSE
1035	CALL ipslerr(3, "grid_toij_scal", "Projection not initilized"," "," ")
1036	ENDIF
1037	!
1038	END SUBROUTINE grid_toij_scal
1039	!
1040	!=========================================================================================
1041	!
1042	SUBROUTINE grid_toij_1d (lon, lat, ri, rj)
1043	!
1044	!
1045	! Argument
1046	REAL(r_std), INTENT(in), DIMENSION(:) :: lon, lat
1047	REAL(r_std), INTENT(out), DIMENSION(:) :: ri, rj
1048	!
1049	! Local
1050	INTEGER :: i, imax
1051	!
1052	imax=SIZE(lon)
1053	!
1054	IF ( proj_stack(1)%code < undef_int ) THEN
1055	DO i=1,imax
1056	!
1057	CALL latlon_to_ij(proj_stack(1), lat(i), lon(i), ri(i), rj(i))
1058	!
1059	ENDDO
1060	ELSE
1061	CALL ipslerr(3, "grid_toij_1d", "Projection not initilized"," "," ")
1062	ENDIF
1063	!
1064	END SUBROUTINE grid_toij_1d
1065	!
1066	!=========================================================================================
1067	!
1068	SUBROUTINE grid_toij_2d (lon, lat, ri, rj)
1069	!
1070	!
1071	! Argument
1072	REAL(r_std), INTENT(in), DIMENSION(:,:) :: lon, lat
1073	REAL(r_std), INTENT(out), DIMENSION(:,:) :: ri, rj
1074	!
1075	! Local
1076	INTEGER :: i, imax, j, jmax
1077	!
1078	imax=SIZE(lon,DIM=1)
1079	jmax=SIZE(lon,DIM=2)
1080	!
1081	IF ( proj_stack(1)%code < undef_int ) THEN
1082	DO i=1,imax
1083	DO j=1,jmax
1084	!
1085	CALL latlon_to_ij(proj_stack(1), lat(i,j), lon(i,j), ri(i,j), rj(i,j))
1086	!
1087	ENDDO
1088	ENDDO
1089	ELSE
1090	CALL ipslerr(3, "grid_toij_2d", "Projection not initilized"," "," ")
1091	ENDIF
1092	!
1093	END SUBROUTINE grid_toij_2d
1094	!
1095	!
1096	!=========================================================================================
1097	!
1098	!
1099	END MODULE grid

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

Download in other formats: