Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

grid.f90 @ 7325

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

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

Context Navigation

source: branches/publications/ORCHIDEE-MICT-BIOENERGY_r7298/src_global/grid.f90 @ 7325

Download in other formats: