Context Navigation

source: branches/publications/ORCHIDEE_CN_CAN_r5698/src_global/grid.f90 @ 7346

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

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

Download in other formats: