Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

grid.f90 @ 8798

Last change on this file since 8798 was 7366, checked in by simon.bowring, 3 years ago
Biochar version
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 36.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	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(out) :: 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, global, &
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, global, &
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	!
659	!
660	IF ( printlev >=4 ) THEN
661	WRITE(numout,*) 'grid_scatter > seglength = ', seglength(1,:)
662	WRITE(numout,*) 'grid_scatter > neighbours = ', neighbours(1,:)
663	WRITE(numout,*) 'grid_scatter > contfrac = ', contfrac(1)
664	WRITE(numout,*) 'grid_scatter > area = ', area(1)
665	ENDIF
666	!
667	END SUBROUTINE grid_scatter
668	!!
669	!!
670	!! =============================================================================================================================
671	!! SUBROUTINE: grid_initproj
672	!!
673	!>\BRIEF Routine to initialise the projection
674	!!
675	!! DESCRIPTION:
676	!!
677	!!
678	!! This subroutine is called by the routine whichs ets-up th grid on which ORCHIDEE is to run.
679	!! The aim is to set-upu the projection so that all the grid variables needed by ORCHIDEE can
680	!! be computed in grid_stuff_regxy
681	!!
682	!! \n
683	!_ ==============================================================================================================================
684	!!
685	!!
686	SUBROUTINE grid_initproj (fid, iim, jjm)
687	!
688	!
689	! 0 interface
690	!
691	IMPLICIT NONE
692	!
693	! 0.1 input !
694	!
695	! Domain size
696	INTEGER(i_std), INTENT(in) :: fid
697	INTEGER(i_std), INTENT(in) :: iim, jjm
698	!
699	! 0.2 Local variables
700	!
701	INTEGER(i_std) :: current_proj, idom, iret, lonid, latid, numLons, numLats
702	INTEGER, DIMENSION(nf90_max_var_dims) :: dimIDs
703	REAL(r_std) :: user_stand_lon, user_truelat1, user_truelat2, user_dxkm, user_dykm
704	REAL(r_std) :: user_dlat, user_dlon, user_known_x, user_known_y, user_known_lat, user_known_lon
705	REAL(r_std), DIMENSION(16) :: corner_lons, corner_lats
706	!
707	INTEGER(i_std) :: iv, i, j
708	CHARACTER(LEN=20) :: varname
709	REAL(r_std) :: dx, dy, dtx, dty, coslat
710	REAL(r_std), ALLOCATABLE, DIMENSION (:) :: LON, LAT
711	REAL(r_std), ALLOCATABLE, DIMENSION (:,:) :: mapfac_x, mapfac_y
712	!
713	!
714	! Only one domain is possible for the moment
715	!
716	idom=1
717	CALL map_init(proj_stack(idom))
718	!
719	! Does ORCHIDEE have the same Earth Radius as the map projection ?
720	!
721	IF ( ABS(R_Earth-EARTH_RADIUS_M) > 0.1 ) THEN
722	WRITE(,) "Earth Radius in WRF : ", EARTH_RADIUS_M
723	WRITE(,) "Earth Radius in ORCHIDEE : ", R_Earth
724	CALL ipslerr (3,'grid_initproj','The Earth radius is not the same in the projection module and ORCHIDEE',&
725	& " ", " ")
726	ENDIF
727	!
728	! Get parameters of the projection from the netCDF file
729	!
730	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "MAP_PROJ", current_proj)
731	!
732	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "STAND_LON", user_stand_lon)
733	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "TRUELAT1", user_truelat1)
734	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "TRUELAT2", user_truelat2)
735	!
736	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "DX", user_dxkm)
737	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "DY", user_dykm)
738	user_dlat = undef
739	user_dlon = undef
740	!
741	IF ( current_proj == PROJ_LATLON ) THEN
742	!
743	iret = NF90_inq_VARID(fid, "XLONG_M",lonid)
744	iret = NF90_INQUIRE_VARIABLE(fid, lonid, dimids = dimIDs)
745	iret = NF90_INQUIRE_DIMENSION(fid, dimIDs(1), len = numLons)
746	iret = NF90_INQUIRE_DIMENSION(fid, dimIDs(2), len = numLats)
747	ALLOCATE(LON(numLons))
748	iret = NF90_GET_VAR(fid, lonid, LON(:), start = (/ 1, 1, 1 /), count = (/ numLons, 1, 1 /))
749
750	iret = NF90_inq_VARID(fid, "XLAT_M",latid)
751	ALLOCATE(LAT(numLats))
752	iret = NF90_GET_VAR(fid, latid, LAT(:), start = (/ 1, 1, 1 /), count = (/ 1, numLats, 1 /))
753
754	user_dlon = (LON(numLons) - LON(1)) / (numLons - 1)
755	user_dlat = (LAT(numLats) - LAT(1)) / (numLats - 1)
756
757	DEALLOCATE(LON,LAT)
758
759	ENDIF
760	! Unable to know from where to get the information
761	user_known_x = 1
762	user_known_y = 1
763	!
764	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "corner_lats", corner_lats)
765	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "corner_lons", corner_lons)
766	user_known_lat = corner_lats(1)
767	user_known_lon = corner_lons(1)
768	!
769	! Read mapfactor, land mask and orography
770	!
771	!
772	! Allocation
773	!
774	ALLOCATE(mapfac_x(iim,jjm))
775	ALLOCATE(mapfac_y(iim,jjm))
776	ALLOCATE(dxwrf(iim,jjm))
777	ALLOCATE(dywrf(iim,jjm))
778	!
779	varname = "MAPFAC_MX"
780	iret = NF90_INQ_VARID (fid, varname, iv)
781	IF (iret /= NF90_NOERR) THEN
782	CALL ipslerr (3,'WRFdomain_Read',"Could not find variable ", varname," ")
783	ELSE
784	iret = NF90_GET_VAR (fid,iv,mapfac_x)
785	ENDIF
786	varname = "MAPFAC_MY"
787	iret = NF90_INQ_VARID (fid, varname, iv)
788	IF (iret /= NF90_NOERR) THEN
789	CALL ipslerr (3,'WRFdomain_Read',"Could not find variable ", varname," ")
790	ELSE
791	iret = NF90_GET_VAR (fid,iv,mapfac_y)
792	ENDIF
793	!
794	! Initilize the projection
795	!
796	if (current_proj == PROJ_LATLON) then
797	call map_set(current_proj, proj_stack(idom), &
798	lat1=user_known_lat, &
799	lon1=user_known_lon, &
800	knowni=user_known_x, &
801	knownj=user_known_y, &
802	latinc=user_dlat, &
803	loninc=user_dlon, &
804	r_earth=R_Earth)
805
806	else if (current_proj == PROJ_MERC) then
807	call map_set(current_proj, proj_stack(idom), &
808	truelat1=user_truelat1, &
809	lat1=user_known_lat, &
810	lon1=user_known_lon, &
811	knowni=user_known_x, &
812	knownj=user_known_y, &
813	dx=user_dxkm, &
814	r_earth=R_Earth)
815
816	else if (current_proj == PROJ_CYL) then
817	call ipslerr(3,"grid_initproj",'Should not have PROJ_CYL as projection for',&
818	'source data in push_source_projection()', " ")
819
820	else if (current_proj == PROJ_CASSINI) then
821	call ipslerr(3,"grid_initproj",'Should not have PROJ_CASSINI as projection for', &
822	'source data in push_source_projection()', " ")
823
824	else if (current_proj == PROJ_LC) then
825	call map_set(current_proj, proj_stack(idom), &
826	truelat1=user_truelat1, &
827	truelat2=user_truelat2, &
828	stdlon=user_stand_lon, &
829	lat1=user_known_lat, &
830	lon1=user_known_lon, &
831	knowni=user_known_x, &
832	knownj=user_known_y, &
833	dx=user_dxkm, &
834	r_earth=R_Earth)
835
836	else if (current_proj == PROJ_ALBERS_NAD83) then
837	call map_set(current_proj, proj_stack(idom), &
838	truelat1=user_truelat1, &
839	truelat2=user_truelat2, &
840	stdlon=user_stand_lon, &
841	lat1=user_known_lat, &
842	lon1=user_known_lon, &
843	knowni=user_known_x, &
844	knownj=user_known_y, &
845	dx=user_dxkm, &
846	r_earth=R_Earth)
847
848	else if (current_proj == PROJ_PS) then
849	call map_set(current_proj, proj_stack(idom), &
850	truelat1=user_truelat1, &
851	stdlon=user_stand_lon, &
852	lat1=user_known_lat, &
853	lon1=user_known_lon, &
854	knowni=user_known_x, &
855	knownj=user_known_y, &
856	dx=user_dxkm, &
857	r_earth=R_Earth)
858
859	else if (current_proj == PROJ_PS_WGS84) then
860	call map_set(current_proj, proj_stack(idom), &
861	truelat1=user_truelat1, &
862	stdlon=user_stand_lon, &
863	lat1=user_known_lat, &
864	lon1=user_known_lon, &
865	knowni=user_known_x, &
866	knownj=user_known_y, &
867	dx=user_dxkm, &
868	r_earth=R_Earth)
869
870	else if (current_proj == PROJ_GAUSS) then
871	call map_set(current_proj, proj_stack(idom), &
872	lat1=user_known_lat, &
873	lon1=user_known_lon, &
874	nlat=nint(user_dlat), &
875	loninc=user_dlon, &
876	r_earth=R_Earth)
877
878	else if (current_proj == PROJ_ROTLL) then
879	call ipslerr(3 ,"grid_initproj",'Should not have PROJ_ROTLL as projection for', &
880	'source data in push_source_projection() as not yet implemented', '')
881	end if
882	!
883	! Transform the mapfactors into dx and dy to be used for the description of the polygons and
884	! interpolations.
885	!
886	DO i=1,iim
887	DO j=1,jjm
888	!
889	IF (proj_stack(idom)%code /= PROJ_LATLON ) THEN
890	dx = proj_stack(idom)%dx
891	! Some projections in WRF do not store dy, in that case dy=dx.
892	IF ( proj_stack(idom)%dy > 0 ) THEN
893	dy = proj_stack(idom)%dy
894	ELSE
895	dy = proj_stack(idom)%dx
896	ENDIF
897	dxwrf(i,j) = dx/mapfac_x(i,j)
898	dywrf(i,j) = dy/mapfac_y(i,j)
899	ELSE
900	!
901	! The LatLon projection is also a special case as here it is not the dx and dy
902	! which are stored in the projection file but the increments in Lon and Lat.
903	!
904	dtx = proj_stack(idom)%loninc
905	dty = proj_stack(idom)%latinc
906	coslat = COS(lat(j) * pi/180. )
907	dxwrf(i,j) = dtx * pi/180. * R_Earth * coslat
908	dywrf(i,j) = dty * pi/180. * R_Earth
909	!
910	ENDIF
911	!
912	ENDDO
913	ENDDO
914	!
915	END SUBROUTINE grid_initproj
916	!
917	!
918	!
919	!=========================================================================================
920	!
921	SUBROUTINE grid_tolola_scal (ri, rj, lon, lat)
922	!
923	!
924	! Argument
925	REAL(r_std), INTENT(in) :: ri, rj
926	REAL(r_std), INTENT(out) :: lon, lat
927	!
928	!
929	IF ( proj_stack(1)%code < undef_int ) THEN
930	!
931	CALL ij_to_latlon(proj_stack(1), ri, rj, lat, lon)
932	!
933	ELSE
934	CALL ipslerr(3, "grid_tolola_scal", "Projection not initilized"," "," ")
935	ENDIF
936	!
937	END SUBROUTINE grid_tolola_scal
938	!
939	!=========================================================================================
940	!
941	SUBROUTINE grid_tolola_1d (ri, rj, lon, lat)
942	!
943	!
944	! Argument
945	REAL(r_std), INTENT(in), DIMENSION(:) :: ri, rj
946	REAL(r_std), INTENT(out), DIMENSION(:) :: lon, lat
947	!
948	! Local
949	INTEGER :: i, imax
950	!
951	imax=SIZE(lon)
952	!
953	IF ( proj_stack(1)%code < undef_int ) THEN
954	DO i=1,imax
955	!
956	CALL ij_to_latlon(proj_stack(1), ri(i), rj(i), lat(i), lon(i))
957	!
958	ENDDO
959	ELSE
960	CALL ipslerr(3, "grid_tolola_1d", "Projection not initilized"," "," ")
961	ENDIF
962	!
963	END SUBROUTINE grid_tolola_1d
964	!
965	!=========================================================================================
966	!
967	SUBROUTINE grid_tolola_2d (ri, rj, lon, lat)
968	!
969	!
970	! Argument
971	REAL(r_std), INTENT(in), DIMENSION(:,:) :: ri, rj
972	REAL(r_std), INTENT(out), DIMENSION(:,:) :: lon, lat
973	!
974	! Local
975	INTEGER :: i, imax, j, jmax
976	!
977	imax=SIZE(lon,DIM=1)
978	jmax=SIZE(lon,DIM=2)
979	!
980	IF ( proj_stack(1)%code < undef_int ) THEN
981	DO i=1,imax
982	DO j=1,jmax
983	!
984	CALL ij_to_latlon(proj_stack(1), ri(i,j), rj(i,j), lat(i,j), lon(i,j))
985	!
986	ENDDO
987	ENDDO
988	ELSE
989	CALL ipslerr(3, "grid_tolola_2d", "Projection not initilized"," "," ")
990	ENDIF
991	!
992	END SUBROUTINE grid_tolola_2d
993	!
994	!=========================================================================================
995	!
996	SUBROUTINE grid_toij_scal (lon, lat, ri, rj)
997	!
998	!
999	! Argument
1000	REAL(r_std), INTENT(in) :: lon, lat
1001	REAL(r_std), INTENT(out) :: ri, rj
1002	!
1003	!
1004	IF ( proj_stack(1)%code < undef_int ) THEN
1005	!
1006	CALL latlon_to_ij(proj_stack(1), lat, lon, ri, rj)
1007	!
1008	ELSE
1009	CALL ipslerr(3, "grid_toij_scal", "Projection not initilized"," "," ")
1010	ENDIF
1011	!
1012	END SUBROUTINE grid_toij_scal
1013	!
1014	!=========================================================================================
1015	!
1016	SUBROUTINE grid_toij_1d (lon, lat, ri, rj)
1017	!
1018	!
1019	! Argument
1020	REAL(r_std), INTENT(in), DIMENSION(:) :: lon, lat
1021	REAL(r_std), INTENT(out), DIMENSION(:) :: ri, rj
1022	!
1023	! Local
1024	INTEGER :: i, imax
1025	!
1026	imax=SIZE(lon)
1027	!
1028	IF ( proj_stack(1)%code < undef_int ) THEN
1029	DO i=1,imax
1030	!
1031	CALL latlon_to_ij(proj_stack(1), lat(i), lon(i), ri(i), rj(i))
1032	!
1033	ENDDO
1034	ELSE
1035	CALL ipslerr(3, "grid_toij_1d", "Projection not initilized"," "," ")
1036	ENDIF
1037	!
1038	END SUBROUTINE grid_toij_1d
1039	!
1040	!=========================================================================================
1041	!
1042	SUBROUTINE grid_toij_2d (lon, lat, ri, rj)
1043	!
1044	!
1045	! Argument
1046	REAL(r_std), INTENT(in), DIMENSION(:,:) :: lon, lat
1047	REAL(r_std), INTENT(out), DIMENSION(:,:) :: ri, rj
1048	!
1049	! Local
1050	INTEGER :: i, imax, j, jmax
1051	!
1052	imax=SIZE(lon,DIM=1)
1053	jmax=SIZE(lon,DIM=2)
1054	!
1055	IF ( proj_stack(1)%code < undef_int ) THEN
1056	DO i=1,imax
1057	DO j=1,jmax
1058	!
1059	CALL latlon_to_ij(proj_stack(1), lat(i,j), lon(i,j), ri(i,j), rj(i,j))
1060	!
1061	ENDDO
1062	ENDDO
1063	ELSE
1064	CALL ipslerr(3, "grid_toij_2d", "Projection not initilized"," "," ")
1065	ENDIF
1066	!
1067	END SUBROUTINE grid_toij_2d
1068	!
1069	!
1070	!=========================================================================================
1071	!
1072	!
1073	END MODULE grid

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

Context Navigation

source: branches/publications/ORCHIDEE_Biochar/src_global/grid.f90 @ 8798

Download in other formats: