Context Navigation

source: branches/publications/ORCHIDEE_DFv1.0_site/src_global/grid.f90 @ 6715

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

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

Download in other formats: