Context Navigation

source: branches/publications/ORCHIDEE_gmd_mict_peat_ch4/src_global/grid.f90 @ 7444

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

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

Download in other formats: