Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

grid.f90 @ 6892

Last change on this file since 6892 was 4977, checked in by simon.bowring, 7 years ago
Currently running (13/02/2018) version includes all necessarily changes to include DOC in MICT code... further parametrisation necessary to equate soil pools with those of normal forcesoil restarts
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 37.2 KB

Line
1
2	!! This module define variables for the grid to gathered points.
3	!!
4	!! @call sechiba_main
5	!! @Version : $Revision$, $Date$
6	!!
7	!< $HeadURL$
8	!< $Date$
9	!< $Revision$
10	!!
11	!! @author Marie-Alice Foujols, Jan Polcher and Martial Mancip
12	!!
13	!! This module archives and makes available for all ORCHIDEE routine the information on the grid
14	!! being used. 3 types of grids are foreseen :
15	!! - Regular longitude latitude grid : This is the default and mostly used for global applications.
16	!! - Regular X/Y grid : this is a typical grid for regional models and requires a projection method
17	!! to go from X/y to lon/lat.
18	!! - unstructures grid : This is a general grid where each cell is a polygone. It prepares ORCHIDEE
19	!! for DYNAMICO.
20	!!
21	!! The subroutines have the following role :
22	!! grid_init : this routine will provide the dimensions needed to allocate the memory and the
23	!! characteristics of the grid.
24	!!
25	!! grid_stuff : This subroutine provides the grid details for all land points. Obviously depending
26	!! on the grid type different level of information need to be provided.
27	!!
28	!f90doc MODULEgrid
29	MODULE grid
30
31	USE grid_var
32	USE defprec
33	USE constantes
34	USE mod_orchidee_para
35
36	USE haversine
37	USE module_llxy
38
39	USE ioipsl
40	USE netcdf
41
42	IMPLICIT NONE
43	!
44	!=================================================================================
45	!
46	! Horizontal grid information
47	!
48	!=================================================================================
49
50	! Global map or not.
51	! There is little chance that if iim <=2 and jjm <= 2 that we have global grid.
52	! Furthermore using the second line allows to avoid pole problems for global grids
53	LOGICAL, SAVE :: global = .TRUE.
54	!$OMP THREADPRIVATE(global)
55
56	! PARAMETERS
57	! default resolution (m)
58	REAL(r_std), PARAMETER :: default_resolution = 250000.
59	!
60	! VARIABLES
61	!
62	!-
63	!- Variable to help describe the grid
64	!- once the points are gathered.
65	!-
66	!! Limits of the domain
67	REAL(r_std), SAVE :: limit_west, limit_east, &
68	& limit_north, limit_south
69	!$OMP THREADPRIVATE(limit_west, limit_east, limit_north, limit_south)
70	!-
71	!! Geographical coordinates
72	REAL(r_std), ALLOCATABLE, DIMENSION (:,:), SAVE :: lalo
73	!$OMP THREADPRIVATE(lalo)
74	!! index of land points
75	INTEGER, ALLOCATABLE, DIMENSION (:), SAVE :: ilandindex,jlandindex
76	!$OMP THREADPRIVATE(ilandindex, jlandindex)
77	!-
78	!! Fraction of continents.
79	REAL(r_std), ALLOCATABLE, DIMENSION (:), SAVE :: contfrac
80	!$OMP THREADPRIVATE(contfrac)
81	!
82	! indices of the NbNeighb neighbours of each grid point
83	! (1=Northern most vertex and then in clockwise order)
84	! Zero or negative index means that this neighbour is not a land point
85	INTEGER(i_std), ALLOCATABLE, DIMENSION (:,:), SAVE :: neighbours
86	!$OMP THREADPRIVATE(neighbours)
87	!
88	! Heading of the direction out of the grid box either through the vertex
89	! of the mid-segment of the polygon.
90	!
91	REAL(r_std), ALLOCATABLE, DIMENSION(:,:), SAVE :: headings
92	!$OMP THREADPRIVATE(headings)
93	!
94	! Length of segments of the polygon.
95	!
96	REAL(r_std), ALLOCATABLE, DIMENSION(:,:), SAVE :: seglength
97	!$OMP THREADPRIVATE(seglength)
98	!
99	! Area of the grid box
100	!
101	REAL(r_std), ALLOCATABLE, DIMENSION(:), SAVE :: area
102	!$OMP THREADPRIVATE(area)
103	!
104	! Coordinats of the vertices
105	!
106	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: corners
107	!$OMP THREADPRIVATE(corners)
108	!
109	! Resolution remains a temporary variable until the merge of the
110	! re-interfacing of the interpolation by Lluis. One this is done
111	! Resolution will be replaced in the model either by area or seglength.
112	!
113	REAL(r_std), ALLOCATABLE, DIMENSION (:,:), SAVE :: resolution
114	!$OMP THREADPRIVATE(resolution)
115	!
116	!
117	!
118	! Get the direction of the grid
119	!
120	CHARACTER(LEN=2), DIMENSION(2), SAVE, PRIVATE :: grid_dir
121	!$OMP THREADPRIVATE(grid_dir)
122	!
123	INTEGER(i_std), PARAMETER :: MAX_DOMAINS=1
124	!
125	type (proj_info), SAVE, dimension(1:MAX_DOMAINS) :: proj_stack
126	!
127	real(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: dxwrf, dywrf
128	!
129	!
130	!=================================================================================
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	CALL bcast(neighbours_g)
450	CALL bcast(resolution_g)
451	!
452	IF ( printlev >= 3 ) WRITE(numout,*) 'Leaving grid_stuff'
453
454	END SUBROUTINE grid_stuff
455	!!
456	!! =============================================================================================================================
457	!! SUBROUTINE: grid_topolylist
458	!!
459	!>\BRIEF This routine transforms a regular grid into a list of polygons which are defined by the following
460	!! quantities :
461	!!
462	!! corners : the n vertices of the polugon in longitude and latitude
463	!! neighbours : the neighbouring land grid box for each of the vertices and segments
464	!! headings : the direction in which the neighbour is
465	!! seglength : the lenght of each segment
466	!! area : the area of the polygon
467	!! ilindex, jlindex : provides the i,j coordinates of the mesh in the global grid.
468	!!
469	!! DESCRIPTION:
470	!!
471	!! \n
472	!_ ==============================================================================================================================
473	!!
474	SUBROUTINE grid_topolylist(gtype, nbseg, nland, iim, jjm, grid_lon, grid_lat, kindex, &
475	& globalg, corners_loc, neighbours_loc, headings_loc, seglength_loc, &
476	& area_loc, ilindex_loc, jlindex_loc)
477	!
478	! 0 interface
479	!
480	IMPLICIT NONE
481	!
482	! 0.1 input !
483	! Grid type
484	CHARACTER(LEN=20), INTENT(in) :: gtype
485	! Number of segments for each polygon
486	INTEGER(i_std), INTENT(in) :: nbseg
487	! Number of land points on the grid
488	INTEGER(i_std), INTENT(in) :: nland
489	! Size of cartesian grid
490	INTEGER(i_std), INTENT(in) :: iim, jjm
491	! Longitudes on cartesian grid
492	REAL(r_std), DIMENSION(iim,jjm), INTENT(in) :: grid_lon
493	! Latitudes on cartesian grid
494	REAL(r_std), DIMENSION(iim,jjm), INTENT(in) :: grid_lat
495	! Index of land point on 2D map (in local position)
496	INTEGER(i_std), DIMENSION(nland), INTENT(in) :: kindex
497	!
498	! 0.2 Output
499	!
500	LOGICAL, INTENT(out) :: globalg
501	!
502	REAL(r_std), DIMENSION(nland,nbseg,2), INTENT(out) :: corners_loc
503	INTEGER(i_std), DIMENSION(nland,nbseg*2), INTENT(out) :: neighbours_loc
504	REAL(r_std), DIMENSION(nland,nbseg*2), INTENT(out) :: headings_loc
505	REAL(r_std), DIMENSION(nland,nbseg), INTENT(out) :: seglength_loc
506	REAL(r_std), DIMENSION(nland), INTENT(out) :: area_loc
507	INTEGER(i_std), DIMENSION(nland), INTENT(out) :: ilindex_loc, jlindex_loc
508	!
509	! 0.3 Local variables
510	!
511	INTEGER(i_std) :: i, is, iss
512	REAL(r_std), DIMENSION(nland,2) :: center
513	REAL(r_std) :: maxdellon, mindellon, maxlon, minlon
514	REAL(r_std), DIMENSION(nland,nbseg*2) :: lonpoly, latpoly
515	!
516	IF ( INDEX(gtype,"RegLonLat") > 0 ) THEN
517	!
518	! If we are in regular Lon Lat, then we test just the longitude and see if we span 0-360deg.
519	!
520	maxdellon=MAXVAL(ABS(grid_lon(1:iim-1,1)-grid_lon(2:iim,1)))
521	mindellon=MINVAL(ABS(grid_lon(1:iim-1,1)-grid_lon(2:iim,1)))
522	maxlon=MAXVAL(grid_lon(1:iim,1))
523	minlon=MINVAL(grid_lon(1:iim,1))
524	!
525	! test if it could be a global grid on 0 -> 360
526	!
527	IF ( minlon > 0 .AND. maxlon > 180 ) THEN
528	IF ( (minlon - maxdellon/2.0 ) <= 0 .AND. (maxlon + maxdellon/2.0) >= 360) THEN
529	globalg = .TRUE.
530	ELSE
531	globalg = .FALSE.
532	ENDIF
533	!
534	! Test if it could be a -180 to 180 grid
535	!
536	ELSE IF ( minlon < 0 .AND. maxlon > 0 ) THEN
537	IF ( (minlon - maxdellon/2.0 ) <= -180 .AND. (maxlon + maxdellon/2.0) >= 180) THEN
538	globalg = .TRUE.
539	ELSE
540	globalg = .FALSE.
541	ENDIF
542	!
543	! If neither condition is met then it cannot be global.
544	!
545	ELSE
546	globalg = .FALSE.
547	ENDIF
548	ELSE IF ( gtype == "RegXY" ) THEN
549	!
550	! The hypothesis is that if we are in RegXY then we are not global
551	!
552	globalg = .FALSE.
553	ELSE
554	STOP "Unknown grid"
555	ENDIF
556	!
557	! 2.0 Transform the grid into a list of polygones while keeping the neighbour relations
558	! between these polygones.
559	!
560	! Each polygone starts with a vertex and alternates vertices and mid-points of segments.
561	!
562	IF (nland == 1) THEN
563	CALL haversine_singlepointploy(iim, jjm, grid_lon, grid_lat, nland, kindex, global, &
564	& nbseg, lonpoly, latpoly, center, &
565	& neighbours_loc, ilindex_loc, jlindex_loc)
566	ELSE IF ( INDEX(gtype, "RegLonLat") > 0 ) THEN
567	CALL haversine_reglatlontoploy(iim, jjm, grid_lon, grid_lat, nland, kindex, global, &
568	& nbseg, lonpoly, latpoly, center, &
569	& neighbours_loc, ilindex_loc, jlindex_loc)
570	ELSE IF ( INDEX(gtype, "RegXY") > 0 ) THEN
571	CALL haversine_regxytoploy(iim, jjm, grid_lon, grid_lat, nland, kindex, proj_stack, &
572	& nbseg, lonpoly, latpoly, center, &
573	& neighbours_loc, ilindex_loc, jlindex_loc)
574	ELSE
575	STOP "Unknown grid"
576	ENDIF
577	!
578	! Save the longitude and latitudes nbseg corners (=vertices) of the polygones
579	!
580	DO i=1,nland
581	DO is=1,nbseg
582	iss=(is-1)*2+1
583	corners_loc(i,is,1) = lonpoly(i,iss)
584	corners_loc(i,is,2) = latpoly(i,iss)
585	ENDDO
586	ENDDO
587	!
588	! Get the heading normal to the 4 segments and through the 4 corners.
589	!
590	CALL haversine_polyheadings(nland, nbseg, lonpoly, latpoly, center, headings_loc)
591	!
592	! Order the points of the polygone in clockwise order Starting with the northern most
593	!
594	CALL haversine_polysort(nland, nbseg, lonpoly, latpoly, headings_loc, neighbours_loc)
595	!
596	! Compute the segment length and area.
597	! For the RegLonLat we have specific calculations for seglength and area.
598	! For projected regular grids we use the great cicle assumption for the segments
599	! but the projected area.
600	! For unstructured grid we use the most general routines.
601	!
602	IF ( INDEX(gtype, "RegLonLat") > 0 ) THEN
603	CALL haversine_laloseglen(nland, nbseg, lonpoly, latpoly, seglength_loc)
604	CALL haversine_laloarea(nland, nbseg, seglength_loc, area_loc)
605	ELSE IF ( INDEX(gtype, "RegXY") > 0 ) THEN
606	CALL haversine_polyseglen(nland, nbseg, lonpoly, latpoly, seglength_loc)
607	CALL haversine_xyarea(nland, nbseg, ilindex_loc, jlindex_loc, dxwrf, dywrf, area_loc)
608	ELSE
609	CALL haversine_polyseglen(nland, nbseg, lonpoly, latpoly, seglength_loc)
610	CALL haversine_polyarea(nland, nbseg, lonpoly, latpoly, area_loc)
611	ENDIF
612	! Compute the area
613
614	!
615	END SUBROUTINE grid_topolylist
616	!!
617	!!
618	!!
619	!! =============================================================================================================================
620	!! SUBROUTINE: grid_scatter
621	!!
622	!>\BRIEF Scatter the grid information so that each processor knows the characteristics of the grid it works on.
623	!!
624	!! DESCRIPTION:
625	!!
626	!!
627	!! The grid information has been computed for the entire grid on the root processor. Now we give each processor
628	!! the information of the piece of the grid it works on. This concerns the following variables describing the grid :
629	!! - area
630	!! - resolution
631	!! - neighbours
632	!! - contfrac : fraction of continent
633	!!
634	!! Should ilandindex and jlandindex not b initialized, we catch-up here. This field is the same on all processors.
635	!!
636	!! TODO :
637	!! This code should get the grid describing fields as arguments and then writem into the *_g variables on
638	!! root_prc before scattering. This would allow to compute the grid characteristics in any subroutine
639	!! fore calling grid_scatter.
640	!!
641	!!
642	!!
643	!! \n
644	!_ ==============================================================================================================================
645	!!
646	!!
647	SUBROUTINE grid_scatter()
648	!
649	!
650	INTEGER(i_std) :: i, ip, jp
651	!
652	IF ( MAXVAL(ilandindex) < 0 .AND. MAXVAL(jlandindex) < 0 ) THEN
653	DO i = 1, nbp_glo
654	!
655	! 1 find numbers of the latitude and longitude of each point
656	!
657
658	! index of latitude
659	jp = INT( (index_g(i)-1) /iim_g ) + 1
660
661	! index of longitude
662	ip = index_g(i) - ( jp-1 ) * iim_g
663	!
664	! Save this information for usage in other modules.
665	!
666	ilandindex(i)=ip
667	jlandindex(i)=jp
668	!
669	ENDDO
670	ENDIF
671	!
672	CALL scatter(neighbours_g, neighbours)
673	CALL scatter(contfrac_g, contfrac)
674	CALL scatter(headings_g, headings)
675	CALL scatter(seglength_g, seglength)
676	CALL scatter(corners_g, corners)
677	CALL scatter(area_g, area)
678	!
679	! TEMPORARY section for resolution
680	!
681	IF ( is_root_prc) THEN
682	IF ( INDEX(GridType,"Reg") > 0 ) THEN
683	resolution_g(:,1) = (seglength_g(:,1)+seglength_g(:,3))/2.0
684	resolution_g(:,2) = (seglength_g(:,2)+seglength_g(:,4))/2.0
685	ELSE
686	CALL ipslerr(3, "grid_scatter", "unsupported grid type.",&
687	& "As long as resolution has not been replaced,",&
688	& "ORCHIDEE cannot run on anything other than regular grids.")
689	ENDIF
690	ENDIF
691	CALL scatter(resolution_g, resolution)
692
693	!
694	!
695	IF ( printlev >=4 ) THEN
696	WRITE(numout,*) 'grid_scatter > seglength = ', seglength(1,:)
697	WRITE(numout,*) 'grid_scatter > neighbours = ', neighbours(1,:)
698	WRITE(numout,*) 'grid_scatter > contfrac = ', contfrac(1)
699	WRITE(numout,*) 'grid_scatter > area = ', area(1)
700	ENDIF
701	!
702	END SUBROUTINE grid_scatter
703	!!
704	!!
705	!! =============================================================================================================================
706	!! SUBROUTINE: grid_initproj
707	!!
708	!>\BRIEF Routine to initialise the projection
709	!!
710	!! DESCRIPTION:
711	!!
712	!!
713	!! This subroutine is called by the routine whichs ets-up th grid on which ORCHIDEE is to run.
714	!! The aim is to set-upu the projection so that all the grid variables needed by ORCHIDEE can
715	!! be computed in grid_stuff_regxy
716	!!
717	!! \n
718	!_ ==============================================================================================================================
719	!!
720	!!
721	SUBROUTINE grid_initproj (fid, iim, jjm)
722	!
723	!
724	! 0 interface
725	!
726	IMPLICIT NONE
727	!
728	! 0.1 input !
729	!
730	! Domain size
731	INTEGER(i_std), INTENT(in) :: fid
732	INTEGER(i_std), INTENT(in) :: iim, jjm
733	!
734	! 0.2 Local variables
735	!
736	INTEGER(i_std) :: current_proj, idom, iret, lonid, latid, numLons, numLats
737	INTEGER, DIMENSION(nf90_max_var_dims) :: dimIDs
738	REAL(r_std) :: user_stand_lon, user_truelat1, user_truelat2, user_dxkm, user_dykm
739	REAL(r_std) :: user_dlat, user_dlon, user_known_x, user_known_y, user_known_lat, user_known_lon
740	REAL(r_std), DIMENSION(16) :: corner_lons, corner_lats
741	!
742	INTEGER(i_std) :: iv, i, j
743	CHARACTER(LEN=20) :: varname
744	REAL(r_std) :: dx, dy, dtx, dty, coslat
745	REAL(r_std), ALLOCATABLE, DIMENSION (:) :: LON, LAT
746	REAL(r_std), ALLOCATABLE, DIMENSION (:,:) :: mapfac_x, mapfac_y
747	!
748	!
749	! Only one domain is possible for the moment
750	!
751	idom=1
752	CALL map_init(proj_stack(idom))
753	!
754	! Does ORCHIDEE have the same Earth Radius as the map projection ?
755	!
756	IF ( ABS(R_Earth-EARTH_RADIUS_M) > 0.1 ) THEN
757	WRITE(,) "Earth Radius in WRF : ", EARTH_RADIUS_M
758	WRITE(,) "Earth Radius in ORCHIDEE : ", R_Earth
759	CALL ipslerr (3,'grid_initproj','The Earth radius is not the same in the projection module and ORCHIDEE',&
760	& " ", " ")
761	ENDIF
762	!
763	! Get parameters of the projection from the netCDF file
764	!
765	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "MAP_PROJ", current_proj)
766	!
767	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "STAND_LON", user_stand_lon)
768	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "TRUELAT1", user_truelat1)
769	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "TRUELAT2", user_truelat2)
770	!
771	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "DX", user_dxkm)
772	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "DY", user_dykm)
773	user_dlat = undef
774	user_dlon = undef
775	!
776	IF ( current_proj == PROJ_LATLON ) THEN
777	!
778	iret = NF90_inq_VARID(fid, "XLONG_M",lonid)
779	iret = NF90_INQUIRE_VARIABLE(fid, lonid, dimids = dimIDs)
780	iret = NF90_INQUIRE_DIMENSION(fid, dimIDs(1), len = numLons)
781	iret = NF90_INQUIRE_DIMENSION(fid, dimIDs(2), len = numLats)
782	ALLOCATE(LON(numLons))
783	iret = NF90_GET_VAR(fid, lonid, LON(:), start = (/ 1, 1, 1 /), count = (/ numLons, 1, 1 /))
784
785	iret = NF90_inq_VARID(fid, "XLAT_M",latid)
786	ALLOCATE(LAT(numLats))
787	iret = NF90_GET_VAR(fid, latid, LAT(:), start = (/ 1, 1, 1 /), count = (/ 1, numLats, 1 /))
788
789	user_dlon = (LON(numLons) - LON(1)) / (numLons - 1)
790	user_dlat = (LAT(numLats) - LAT(1)) / (numLats - 1)
791
792	DEALLOCATE(LON,LAT)
793
794	ENDIF
795	! Unable to know from where to get the information
796	user_known_x = 1
797	user_known_y = 1
798	!
799	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "corner_lats", corner_lats)
800	iret = NF90_GET_ATT(fid, NF90_GLOBAL, "corner_lons", corner_lons)
801	user_known_lat = corner_lats(1)
802	user_known_lon = corner_lons(1)
803	!
804	! Read mapfactor, land mask and orography
805	!
806	!
807	! Allocation
808	!
809	ALLOCATE(mapfac_x(iim,jjm))
810	ALLOCATE(mapfac_y(iim,jjm))
811	ALLOCATE(dxwrf(iim,jjm))
812	ALLOCATE(dywrf(iim,jjm))
813	!
814	varname = "MAPFAC_MX"
815	iret = NF90_INQ_VARID (fid, varname, iv)
816	IF (iret /= NF90_NOERR) THEN
817	CALL ipslerr (3,'WRFdomain_Read',"Could not find variable ", varname," ")
818	ELSE
819	iret = NF90_GET_VAR (fid,iv,mapfac_x)
820	ENDIF
821	varname = "MAPFAC_MY"
822	iret = NF90_INQ_VARID (fid, varname, iv)
823	IF (iret /= NF90_NOERR) THEN
824	CALL ipslerr (3,'WRFdomain_Read',"Could not find variable ", varname," ")
825	ELSE
826	iret = NF90_GET_VAR (fid,iv,mapfac_y)
827	ENDIF
828	!
829	! Initilize the projection
830	!
831	if (current_proj == PROJ_LATLON) then
832	call map_set(current_proj, proj_stack(idom), &
833	lat1=user_known_lat, &
834	lon1=user_known_lon, &
835	knowni=user_known_x, &
836	knownj=user_known_y, &
837	latinc=user_dlat, &
838	loninc=user_dlon, &
839	r_earth=R_Earth)
840
841	else if (current_proj == PROJ_MERC) then
842	call map_set(current_proj, proj_stack(idom), &
843	truelat1=user_truelat1, &
844	lat1=user_known_lat, &
845	lon1=user_known_lon, &
846	knowni=user_known_x, &
847	knownj=user_known_y, &
848	dx=user_dxkm, &
849	r_earth=R_Earth)
850
851	else if (current_proj == PROJ_CYL) then
852	call ipslerr(3,"grid_initproj",'Should not have PROJ_CYL as projection for',&
853	'source data in push_source_projection()', " ")
854
855	else if (current_proj == PROJ_CASSINI) then
856	call ipslerr(3,"grid_initproj",'Should not have PROJ_CASSINI as projection for', &
857	'source data in push_source_projection()', " ")
858
859	else if (current_proj == PROJ_LC) then
860	call map_set(current_proj, proj_stack(idom), &
861	truelat1=user_truelat1, &
862	truelat2=user_truelat2, &
863	stdlon=user_stand_lon, &
864	lat1=user_known_lat, &
865	lon1=user_known_lon, &
866	knowni=user_known_x, &
867	knownj=user_known_y, &
868	dx=user_dxkm, &
869	r_earth=R_Earth)
870
871	else if (current_proj == PROJ_ALBERS_NAD83) then
872	call map_set(current_proj, proj_stack(idom), &
873	truelat1=user_truelat1, &
874	truelat2=user_truelat2, &
875	stdlon=user_stand_lon, &
876	lat1=user_known_lat, &
877	lon1=user_known_lon, &
878	knowni=user_known_x, &
879	knownj=user_known_y, &
880	dx=user_dxkm, &
881	r_earth=R_Earth)
882
883	else if (current_proj == PROJ_PS) then
884	call map_set(current_proj, proj_stack(idom), &
885	truelat1=user_truelat1, &
886	stdlon=user_stand_lon, &
887	lat1=user_known_lat, &
888	lon1=user_known_lon, &
889	knowni=user_known_x, &
890	knownj=user_known_y, &
891	dx=user_dxkm, &
892	r_earth=R_Earth)
893
894	else if (current_proj == PROJ_PS_WGS84) then
895	call map_set(current_proj, proj_stack(idom), &
896	truelat1=user_truelat1, &
897	stdlon=user_stand_lon, &
898	lat1=user_known_lat, &
899	lon1=user_known_lon, &
900	knowni=user_known_x, &
901	knownj=user_known_y, &
902	dx=user_dxkm, &
903	r_earth=R_Earth)
904
905	else if (current_proj == PROJ_GAUSS) then
906	call map_set(current_proj, proj_stack(idom), &
907	lat1=user_known_lat, &
908	lon1=user_known_lon, &
909	nlat=nint(user_dlat), &
910	loninc=user_dlon, &
911	r_earth=R_Earth)
912
913	else if (current_proj == PROJ_ROTLL) then
914	call ipslerr(3 ,"grid_initproj",'Should not have PROJ_ROTLL as projection for', &
915	'source data in push_source_projection() as not yet implemented', '')
916	end if
917	!
918	! Transform the mapfactors into dx and dy to be used for the description of the polygons and
919	! interpolations.
920	!
921	DO i=1,iim
922	DO j=1,jjm
923	!
924	IF (proj_stack(idom)%code /= PROJ_LATLON ) THEN
925	dx = proj_stack(idom)%dx
926	! Some projections in WRF do not store dy, in that case dy=dx.
927	IF ( proj_stack(idom)%dy > 0 ) THEN
928	dy = proj_stack(idom)%dy
929	ELSE
930	dy = proj_stack(idom)%dx
931	ENDIF
932	dxwrf(i,j) = dx/mapfac_x(i,j)
933	dywrf(i,j) = dy/mapfac_y(i,j)
934	ELSE
935	!
936	! The LatLon projection is also a special case as here it is not the dx and dy
937	! which are stored in the projection file but the increments in Lon and Lat.
938	!
939	dtx = proj_stack(idom)%loninc
940	dty = proj_stack(idom)%latinc
941	coslat = COS(lat(j) * pi/180. )
942	dxwrf(i,j) = dtx * pi/180. * R_Earth * coslat
943	dywrf(i,j) = dty * pi/180. * R_Earth
944	!
945	ENDIF
946	!
947	ENDDO
948	ENDDO
949	!
950	END SUBROUTINE grid_initproj
951	!
952	!
953	!
954	!=========================================================================================
955	!
956	SUBROUTINE grid_tolola_scal (ri, rj, lon, lat)
957	!
958	!
959	! Argument
960	REAL(r_std), INTENT(in) :: ri, rj
961	REAL(r_std), INTENT(out) :: lon, lat
962	!
963	!
964	IF ( proj_stack(1)%code < undef_int ) THEN
965	!
966	CALL ij_to_latlon(proj_stack(1), ri, rj, lat, lon)
967	!
968	ELSE
969	CALL ipslerr(3, "grid_tolola_scal", "Projection not initilized"," "," ")
970	ENDIF
971	!
972	END SUBROUTINE grid_tolola_scal
973	!
974	!=========================================================================================
975	!
976	SUBROUTINE grid_tolola_1d (ri, rj, lon, lat)
977	!
978	!
979	! Argument
980	REAL(r_std), INTENT(in), DIMENSION(:) :: ri, rj
981	REAL(r_std), INTENT(out), DIMENSION(:) :: lon, lat
982	!
983	! Local
984	INTEGER :: i, imax
985	!
986	imax=SIZE(lon)
987	!
988	IF ( proj_stack(1)%code < undef_int ) THEN
989	DO i=1,imax
990	!
991	CALL ij_to_latlon(proj_stack(1), ri(i), rj(i), lat(i), lon(i))
992	!
993	ENDDO
994	ELSE
995	CALL ipslerr(3, "grid_tolola_1d", "Projection not initilized"," "," ")
996	ENDIF
997	!
998	END SUBROUTINE grid_tolola_1d
999	!
1000	!=========================================================================================
1001	!
1002	SUBROUTINE grid_tolola_2d (ri, rj, lon, lat)
1003	!
1004	!
1005	! Argument
1006	REAL(r_std), INTENT(in), DIMENSION(:,:) :: ri, rj
1007	REAL(r_std), INTENT(out), DIMENSION(:,:) :: lon, lat
1008	!
1009	! Local
1010	INTEGER :: i, imax, j, jmax
1011	!
1012	imax=SIZE(lon,DIM=1)
1013	jmax=SIZE(lon,DIM=2)
1014	!
1015	IF ( proj_stack(1)%code < undef_int ) THEN
1016	DO i=1,imax
1017	DO j=1,jmax
1018	!
1019	CALL ij_to_latlon(proj_stack(1), ri(i,j), rj(i,j), lat(i,j), lon(i,j))
1020	!
1021	ENDDO
1022	ENDDO
1023	ELSE
1024	CALL ipslerr(3, "grid_tolola_2d", "Projection not initilized"," "," ")
1025	ENDIF
1026	!
1027	END SUBROUTINE grid_tolola_2d
1028	!
1029	!=========================================================================================
1030	!
1031	SUBROUTINE grid_toij_scal (lon, lat, ri, rj)
1032	!
1033	!
1034	! Argument
1035	REAL(r_std), INTENT(in) :: lon, lat
1036	REAL(r_std), INTENT(out) :: ri, rj
1037	!
1038	!
1039	IF ( proj_stack(1)%code < undef_int ) THEN
1040	!
1041	CALL latlon_to_ij(proj_stack(1), lat, lon, ri, rj)
1042	!
1043	ELSE
1044	CALL ipslerr(3, "grid_toij_scal", "Projection not initilized"," "," ")
1045	ENDIF
1046	!
1047	END SUBROUTINE grid_toij_scal
1048	!
1049	!=========================================================================================
1050	!
1051	SUBROUTINE grid_toij_1d (lon, lat, ri, rj)
1052	!
1053	!
1054	! Argument
1055	REAL(r_std), INTENT(in), DIMENSION(:) :: lon, lat
1056	REAL(r_std), INTENT(out), DIMENSION(:) :: ri, rj
1057	!
1058	! Local
1059	INTEGER :: i, imax
1060	!
1061	imax=SIZE(lon)
1062	!
1063	IF ( proj_stack(1)%code < undef_int ) THEN
1064	DO i=1,imax
1065	!
1066	CALL latlon_to_ij(proj_stack(1), lat(i), lon(i), ri(i), rj(i))
1067	!
1068	ENDDO
1069	ELSE
1070	CALL ipslerr(3, "grid_toij_1d", "Projection not initilized"," "," ")
1071	ENDIF
1072	!
1073	END SUBROUTINE grid_toij_1d
1074	!
1075	!=========================================================================================
1076	!
1077	SUBROUTINE grid_toij_2d (lon, lat, ri, rj)
1078	!
1079	!
1080	! Argument
1081	REAL(r_std), INTENT(in), DIMENSION(:,:) :: lon, lat
1082	REAL(r_std), INTENT(out), DIMENSION(:,:) :: ri, rj
1083	!
1084	! Local
1085	INTEGER :: i, imax, j, jmax
1086	!
1087	imax=SIZE(lon,DIM=1)
1088	jmax=SIZE(lon,DIM=2)
1089	!
1090	IF ( proj_stack(1)%code < undef_int ) THEN
1091	DO i=1,imax
1092	DO j=1,jmax
1093	!
1094	CALL latlon_to_ij(proj_stack(1), lat(i,j), lon(i,j), ri(i,j), rj(i,j))
1095	!
1096	ENDDO
1097	ENDDO
1098	ELSE
1099	CALL ipslerr(3, "grid_toij_2d", "Projection not initilized"," "," ")
1100	ENDIF
1101	!
1102	END SUBROUTINE grid_toij_2d
1103	!
1104	!
1105	!=========================================================================================
1106	!
1107	!
1108	END MODULE grid

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

Context Navigation

source: branches/publications/ORCHIDEE_gmd_2018_MICT-LEAK/src_global/grid.f90 @ 6892

Download in other formats: