Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

grid.f90 @ 7541

Last change on this file since 7541 was 7541, checked in by fabienne.maignan, 2 years ago
Zhang publication on coupling factor
File size: 38.0 KB

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

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

Context Navigation

source: branches/publications/ORCHIDEE_2.2_r7266/ORCHIDEE/src_global/grid.f90 @ 7541

Download in other formats: