source: branches/publications/ORCHIDEE_gmd-2018-182/src_global/grid.f90 @ 5444

Last change on this file since 5444 was 1078, checked in by anne.cozic, 12 years ago

Merge between branche OpenMP2 at revision 1076 and trunk revision 1062

this merge doesn't change results for Orchidee with compilation MPI

test with OFFLINE and LMDZOR

There is still a bug when the modele is compile with OpenMP
  • Property svn:keywords set to Revision Date HeadURL Date Author Revision
File size: 17.5 KB
Line 
1
2!! This module define variables for the grid to gathered points.
3!!
4!! @call sechiba_main
5!! @Version : $Revision$, $Date$
6!!
7!< $HeadURL$
8!< $Date$
9!< $Author$
10!< $Revision$
11!!
12!! @author Marie-Alice Foujols, Jan Polcher and Martial Mancip
13!!
14!!
15!f90doc MODULEgrid
16MODULE grid
17
18  USE defprec
19  USE constantes
20  USE mod_orchidee_para
21
22  IMPLICIT NONE
23
24  !
25  ! PARAMETERS
26  ! default resolution (m)
27  REAL(r_std), PARAMETER :: default_resolution = 250000.
28  !
29  ! VARIABLES
30  !
31  ! Global map or not.
32  ! There is little change that if iim <=2 and jjm <= 2 that we have global grid.
33  ! Furthermore using the second line allows to avoid pole problems for global grids
34  LOGICAL, SAVE                                      :: global = .TRUE.
35!$OMP THREADPRIVATE(global)
36  !
37  !-
38  !- Variable to help describe the grid
39  !- once the points are gathered.
40  !-
41  !! Limits of the domain
42  REAL(r_std), SAVE                                   :: limit_west, limit_east, &
43       &                                                limit_north, limit_south
44!$OMP THREADPRIVATE(limit_west, limit_east, limit_north, limit_south)
45  !-
46  !! Geographical coordinates
47  REAL(r_std), ALLOCATABLE, DIMENSION (:,:), SAVE     :: lalo
48!$OMP THREADPRIVATE(lalo)
49  !! index  of land points
50  INTEGER, ALLOCATABLE, DIMENSION (:), SAVE          :: ilandindex,jlandindex
51!$OMP THREADPRIVATE(ilandindex, jlandindex)
52  !-
53  !! Fraction of continents.
54  REAL(r_std), ALLOCATABLE, DIMENSION (:), SAVE       :: contfrac
55!$OMP THREADPRIVATE(contfrac)
56  !
57  ! indices of the 4 neighbours of each grid point (1=N, 2=E, 3=S, 4=W)
58  !   a zero or negative index means that this neighbour is not a land point
59  INTEGER(i_std), ALLOCATABLE, DIMENSION (:,:), SAVE :: neighbours
60!$OMP THREADPRIVATE(neighbours)
61  !
62  ! resolution at each grid point in m (1=E-W, 2=N-S)
63  ! (size in x an y of the grid)
64  REAL(r_std), ALLOCATABLE, DIMENSION (:,:), SAVE     :: resolution
65!$OMP THREADPRIVATE(resolution)
66  REAL(r_std), DIMENSION(2), SAVE :: min_resol,max_resol
67!$OMP THREADPRIVATE(min_resol, max_resol)
68  REAL(r_std), ALLOCATABLE, DIMENSION (:), SAVE     :: area
69!$OMP THREADPRIVATE(area)
70  !
71  !
72  ! Get the direction of the grid
73  !
74  CHARACTER(LEN=2), DIMENSION(2), SAVE, PRIVATE      :: grid_dir
75!$OMP THREADPRIVATE(grid_dir)
76  !
77  ! Rose gives the geographical direction for the various index increments
78  ! The following corespondences exist
79  !                       WE&NS  WE&SN and so on !
80  ! rose(1) = i+0 & j-1    NN     SS
81  ! rose(2) = i+1 & j-1    NE     SE
82  ! rose(3) = i+1 & j+0    EE     EE
83  ! rose(4) = i+1 & j+1    SE     NE
84  ! rose(5) = i+0 & j+1    SS     NN
85  ! rose(6) = i-1 & j+1    SW     NW
86  ! rose(7) = i-1 & j+0    WW     WW
87  ! rose(8) = i-1 & j-1    NW     SW
88  INTEGER(i_std), DIMENSION(8), SAVE, PRIVATE        :: rose
89!$OMP THREADPRIVATE(rose)
90  !
91  ! The calendar
92  CHARACTER(LEN=20), SAVE         :: calendar_str
93!$OMP THREADPRIVATE(calendar_str)
94  !
95  ! The date
96  REAL(r_std), SAVE                     :: in_julian
97!$OMP THREADPRIVATE(in_julian)
98  ! Diff with day 0
99  REAL(r_std), SAVE                     :: julian_diff
100!$OMP THREADPRIVATE(julian_diff)
101  !
102  INTEGER(i_std), SAVE                  :: year, month, day
103!$OMP THREADPRIVATE(year, month, day)
104  REAL(r_std), SAVE                     :: sec
105!$OMP THREADPRIVATE(sec)
106  !
107  ! month_len (d)
108  INTEGER(i_std), SAVE                  :: month_len
109!$OMP THREADPRIVATE(month_len)
110  !
111  ! year length (d)
112  INTEGER(i_std), SAVE            :: year_length=0
113!$OMP THREADPRIVATE(year_length)
114  !
115  ! Ration between calendar year in days (ie 360d or 365d ...) to gregorian year length
116  REAL(r_std), SAVE :: year_spread
117!$OMP THREADPRIVATE(year_spread)
118  !
119CONTAINS
120  !
121  !f90doc CONTAINS
122  !
123  !
124  SUBROUTINE init_grid ( npts )
125    !
126    ! 0 interface
127    !
128    IMPLICIT NONE
129    !
130    ! 0.1 input  !
131   
132    ! Domain size
133    INTEGER(i_std), INTENT(in)                                 :: npts !! Number of local continental points
134    !
135    !  Create the internal coordinate table
136    !
137    IF ( (.NOT.ALLOCATED(lalo))) THEN
138       ALLOCATE(lalo(npts,2))
139       lalo(:,:) = val_exp
140    ENDIF
141    !-
142    !- Store variable to help describe the grid
143    !- once the points are gathered.
144    !-
145    IF ( (.NOT.ALLOCATED(neighbours))) THEN
146       ALLOCATE(neighbours(npts,8))
147       neighbours(:,:) = -999999
148    ENDIF
149    IF ( (.NOT.ALLOCATED(resolution))) THEN
150       ALLOCATE(resolution(npts,2))
151       resolution(:,:) = val_exp
152    ENDIF
153    IF ( (.NOT.ALLOCATED(area))) THEN
154       ALLOCATE(area(npts))
155       area(:) = val_exp
156    ENDIF
157    !
158    !- Store the fraction of the continents only once so that the user
159    !- does not change them afterwards.
160    !
161    IF ( (.NOT.ALLOCATED(contfrac))) THEN
162       ALLOCATE(contfrac(npts))
163       contfrac(:) = val_exp
164    ENDIF
165    !
166    ! Allocation of index coordinates
167    IF (.NOT. ALLOCATED(ilandindex)) THEN
168       ALLOCATE(ilandindex(npts),jlandindex(npts))
169       ilandindex(:) = -10000000
170       jlandindex(:) = -10000000
171    ENDIF
172    !
173  END SUBROUTINE init_grid
174
175  SUBROUTINE grid_stuff (npts_glo, iim, jjm, grid_lon, grid_lat, kindex)
176    !
177    ! 0 interface
178    !
179    IMPLICIT NONE
180    !
181    ! 0.1 input  !
182   
183    ! Domain size
184    INTEGER(i_std), INTENT(in)                                 :: npts_glo
185    ! Size of cartesian grid
186    INTEGER(i_std), INTENT(in)                                 :: iim, jjm
187    ! Longitudes on cartesian grid
188    REAL(r_std), DIMENSION(iim,jjm), INTENT(in)                 :: grid_lon
189    ! Latitudes on cartesian grid
190    REAL(r_std), DIMENSION(iim,jjm), INTENT(in)                 :: grid_lat
191    ! Index of land point on 2D map (in local position)
192    INTEGER(i_std), DIMENSION(:), INTENT(in)                :: kindex   
193    !
194    ! 0.3 local
195    !
196    ! Index of land point on 2D map (in global position)
197    INTEGER, ALLOCATABLE, DIMENSION (:)                        :: index_p
198    !
199    ! which STOMATE point corresponds to the given point on the cartesian grid
200    INTEGER(i_std), DIMENSION(iim,jjm)                   :: correspondance
201    ! cosine of the latitude
202    REAL(r_std)                                           :: coslat
203    ! number of points where default resolution is used
204    INTEGER(i_std)                                       :: ndefault_lon, ndefault_lat
205    ! Indices
206    INTEGER(i_std)                                       :: i,ip,jp, imm1, imp1, imm1l, imp1l, ii
207    !
208    INTEGER(i_std), SAVE                           :: bavard=2
209!$OMP THREADPRIVATE(bavard)
210   
211    ! =========================================================================
212   
213    IF ( bavard .GE. 4 ) WRITE(numout,*) 'Entering grid_stuff'
214
215    ! default resolution
216    IF ( bavard .GT. 1 ) WRITE(numout,*) 'grid stuff: default resolution (m): ',default_resolution
217    !
218    !-
219    IF (is_root_prc) THEN
220       ! Check if we have a global map or not.
221       ! There is little change that if iim <=2 and jjm <= 2 that we have global grid.
222       ! Furthermore using the second line allows to avoid pole problems for global grids
223       IF (iim <= 2 .OR. jjm <= 2) THEN
224          global = .FALSE.
225       ELSE
226          ! We assume here that the longitude is in increasing order and in degrees.
227          IF ( grid_lon(iim,2)-grid_lon(1,2) >= 360. - (grid_lon(2,2)-grid_lon(1,2)) ) THEN
228             global = .TRUE.
229          ELSE
230             global = .FALSE.
231          ENDIF
232       ENDIF
233       !
234       ! Get the direction of the grid
235       !
236       IF ( iim > 1 ) THEN
237          IF ( grid_lon(1,1) <= grid_lon(2,1) ) THEN
238             grid_dir(1) = 'WE'
239          ELSE
240             grid_dir(1) = 'EW'
241          ENDIF
242       ELSE
243          grid_dir(1) = 'WE'
244       ENDIF
245       !
246       IF ( jjm > 1 ) THEN
247          IF ( grid_lat(1,1) >= grid_lat(1,2) ) THEN
248             grid_dir(2) = 'NS'
249          ELSE
250             grid_dir(2) = 'SN'
251          ENDIF
252       ELSE
253          grid_dir(2) = 'NS'
254       ENDIF
255       !
256       !! WRITE(numout,*) 'Longitude direction :', grid_dir(1)
257       !! WRITE(numout,*) 'Latitude  direction :', grid_dir(2)
258       !
259       ndefault_lon = 0
260       ndefault_lat = 0
261       ! initialize output
262       neighbours_g(:,:) = -1
263       resolution_g(:,:) = zero
264       min_resol(:) = 1.e6
265       max_resol(:) = moins_un
266       
267       correspondance(:,:) = -1
268       DO i = 1, npts_glo         
269          !
270          ! 1 find numbers of the latitude and longitude of each point
271          !
272         
273          ! index of latitude
274          jp = INT( (index_g(i)-1) /iim ) + 1
275         
276          ! index of longitude
277          ip = index_g(i) - ( jp-1 ) * iim
278          !
279          !correspondance(ip,jp) = kindex(i)
280          !
281          correspondance(ip,jp) = i
282
283       ENDDO
284 
285       !
286       ! Get the "wind rose" for the various orientation of the grid
287       !
288       IF ( grid_dir(1) .EQ. 'WE' .AND.  grid_dir(2) .EQ. 'NS' ) THEN
289          rose(1) = 1
290          rose(2) = 2
291          rose(3) = 3
292          rose(4) = 4
293          rose(5) = 5
294          rose(6) = 6
295          rose(7) = 7
296          rose(8) = 8
297       ELSE IF ( grid_dir(1) .EQ. 'EW' .AND.  grid_dir(2) .EQ. 'NS' ) THEN
298          rose(1) = 1
299          rose(2) = 8
300          rose(3) = 7
301          rose(4) = 6
302          rose(5) = 5
303          rose(6) = 4
304          rose(7) = 3
305          rose(8) = 2
306       ELSE IF ( grid_dir(1) .EQ. 'WE' .AND.  grid_dir(2) .EQ. 'SN' ) THEN
307          rose(1) = 5
308          rose(2) = 4
309          rose(3) = 3
310          rose(4) = 2
311          rose(5) = 1
312          rose(6) = 8
313          rose(7) = 7
314          rose(8) = 6
315       ELSE IF ( grid_dir(1) .EQ. 'EW' .AND.  grid_dir(2) .EQ. 'SN' ) THEN
316          rose(1) = 5
317          rose(2) = 6
318          rose(3) = 7
319          rose(4) = 8
320          rose(5) = 1
321          rose(6) = 2
322          rose(7) = 3
323          rose(8) = 4
324       ELSE
325          WRITE(numout,*) 'We can not be here'
326          STOP 'grid_stuff'
327       ENDIF
328
329       DO i = 1, npts_glo
330
331          ! index of latitude
332          jp = INT( (index_g(i)-1) /iim ) + 1
333         
334          ! index of longitude
335          ip = index_g(i) - ( jp-1 ) * iim
336         
337          !
338          ! 2 resolution
339          !
340         
341          !
342          ! 2.1 longitude
343          !
344         
345          ! prevent infinite resolution at the pole
346          coslat = MAX( COS( grid_lat(ip,jp) * pi/180. ), mincos )     
347          IF ( iim .GT. 1 ) THEN
348             
349             IF ( ip .EQ. 1 ) THEN
350                resolution_g(i,1) = &
351                     ABS( grid_lon(ip+1,jp) - grid_lon(ip,jp) ) * &
352                     pi/180. * R_Earth * coslat
353             ELSEIF ( ip .EQ. iim ) THEN
354                resolution_g(i,1) = &
355                     ABS( grid_lon(ip,jp) - grid_lon(ip-1,jp) ) * &
356                     pi/180. * R_Earth * coslat
357             ELSE
358                resolution_g(i,1) = &
359                     ABS( grid_lon(ip+1,jp) - grid_lon(ip-1,jp) )/2. *&
360                     pi/180. * R_Earth * coslat
361             ENDIF
362             
363          ELSE
364             
365             resolution_g(i,1) = default_resolution
366             
367             ndefault_lon = ndefault_lon + 1
368
369          ENDIF
370
371          !
372          ! 2.2 latitude
373          !
374         
375          IF ( jjm .GT. 1 ) THEN
376
377             IF ( jp .EQ. 1 ) THEN
378                resolution_g(i,2) = &
379                     ABS( grid_lat(ip,jp) - grid_lat(ip,jp+1) ) * &
380                     pi/180. * R_Earth
381             ELSEIF ( jp .EQ. jjm ) THEN
382                resolution_g(i,2) = &
383                     ABS( grid_lat(ip,jp-1) - grid_lat(ip,jp) ) * &
384                     pi/180. * R_Earth
385             ELSE
386                resolution_g(i,2) = &
387                     ABS( grid_lat(ip,jp-1) - grid_lat(ip,jp+1) )/2. *&
388                     pi/180. * R_Earth
389             ENDIF
390             
391          ELSE
392             
393             resolution_g(i,2) = default_resolution
394             
395             ndefault_lat = ndefault_lat + 1
396             
397          ENDIF
398          min_resol(1) = MIN(resolution_g(i,1),min_resol(1))
399          min_resol(2) = MIN(resolution_g(i,2),min_resol(2))
400          max_resol(1) = MAX(resolution_g(i,1),max_resol(1))
401          max_resol(2) = MAX(resolution_g(i,2),max_resol(2))
402
403          area_g(i) = resolution_g(i,1)*resolution_g(i,2)
404
405          !
406          ! 3 find neighbours
407          !     
408          imm1 = 0
409          IF ( ip .GT. 1 ) THEN
410             imm1 = ip - 1
411          ELSEIF ( global ) THEN
412             imm1 = iim
413          ENDIF
414         
415          imp1 = 0
416          IF ( ip .LT. iim ) THEN
417             imp1 = ip + 1
418          ELSEIF ( global ) THEN
419             imp1 = 1
420          ENDIF
421          !
422          ! East and West
423          !
424          IF ( imp1 > 0 ) THEN
425             neighbours_g(i,rose(3)) = correspondance(imp1,jp)
426          ELSE
427             neighbours_g(i,rose(3)) = -1
428          ENDIF
429          IF ( imm1 > 0 ) THEN
430             neighbours_g(i,rose(7)) = correspondance(imm1,jp)
431          ELSE
432             neighbours_g(i,rose(7)) = -1
433          ENDIF
434          !
435          ! North
436          !
437          IF ( jp .GT. 1 ) THEN
438
439             neighbours_g(i,rose(1)) = correspondance(ip,jp-1)
440             IF ( imp1 > 0 ) THEN
441                neighbours_g(i,rose(2)) = correspondance(imp1,jp-1)
442             ELSE
443                neighbours_g(i,rose(2)) = -1
444             ENDIF
445             IF ( imm1 > 0 ) THEN
446                neighbours_g(i,rose(8)) = correspondance(imm1,jp-1)
447             ELSE
448                neighbours_g(i,rose(8)) = -1
449             ENDIF
450
451          ELSE
452             IF ( global ) THEN
453               
454                ! special treatment for the pole if we are really in a 2d grid
455               
456                IF ( ( iim .GT. 1 ) .AND. ( jjm .GT. 1 ) ) THEN
457                   !
458                   ii = MOD(ip+iim/2-1,iim)+1
459                   imm1l = ii - 1
460                   IF ( imm1l .LT. 1 ) imm1l = iim
461                   imp1l = ii + 1
462                   IF ( imp1l .GT. iim ) imp1l = 1
463                   !
464                   IF ( ABS( ( grid_lat(ip,jp) ) - 90. ) .LT. min_sechiba ) THEN
465                      ! the grid point sits exactly on the pole. The neighbour is situated
466                   ! at a lower latitude.
467                      neighbours_g(i,rose(1)) = correspondance( ii, jp+1 )
468                      neighbours_g(i,rose(2)) = correspondance( imm1l, jp+1 )
469                      neighbours_g(i,rose(8)) = correspondance( imp1l, jp+1 )
470                   ELSE
471                      ! look across the North Pole
472                      neighbours_g(i,rose(1)) = correspondance( ii, jp )
473                      neighbours_g(i,rose(2)) = correspondance( imm1l, jp )
474                      neighbours_g(i,rose(8)) = correspondance( imp1l, jp )
475                   ENDIF
476                ENDIF
477
478             ELSE
479
480                neighbours_g(i,rose(1)) = -1
481                neighbours_g(i,rose(2)) = -1
482                neighbours_g(i,rose(8)) = -1
483               
484             ENDIF
485
486          ENDIF
487         
488          ! South
489          IF ( jp .LT. jjm ) THEN
490
491             neighbours_g(i,rose(5)) = correspondance(ip,jp+1)
492             IF ( imp1 > 0 ) THEN
493                neighbours_g(i,rose(4)) = correspondance(imp1,jp+1)
494             ELSE
495                neighbours_g(i,rose(4)) = -1
496             ENDIF
497             IF ( imm1 > 0 ) THEN
498                neighbours_g(i,rose(6)) = correspondance(imm1,jp+1)
499             ELSE
500                neighbours_g(i,rose(6)) = -1
501             ENDIF
502
503          ELSE
504             
505             IF ( global ) THEN
506
507                ! special treatment for the pole if we are really in a 2d grid
508               
509                IF ( ( iim .GT. 1 ) .AND. ( jjm .GT. 1 ) ) THEN
510                   !
511                   ii = MOD(ip+iim/2-1,iim)+1
512                   imm1l = ii - 1
513                   IF ( imm1l .LT. 1 ) imm1l = iim
514                   imp1l = ii + 1
515                   IF ( imp1l .GT. iim ) imp1l = 1
516                   !
517                   IF ( ( ABS( grid_lat(ip,jp) ) - 90. ) .LT. min_sechiba ) THEN
518                      ! the grid point sits exactly on the pole. The neighbour is situated
519                      ! at a lower latitude.
520                      neighbours_g(i,rose(5)) = correspondance( ii, jp-1 )
521                      neighbours_g(i,rose(4)) = correspondance( imm1l, jp-1 )
522                      neighbours_g(i,rose(6)) = correspondance( imp1l, jp-1 )
523                   ELSE
524                      ! look across the South Pole
525                      neighbours_g(i,rose(5)) = correspondance( ii, jp )
526                      neighbours_g(i,rose(4)) = correspondance( imm1l, jp )
527                      neighbours_g(i,rose(6)) = correspondance( imp1l, jp )
528                   ENDIF
529                ENDIF
530
531             ELSE
532               
533                neighbours_g(i,rose(5)) = -1
534                neighbours_g(i,rose(4)) = -1
535                neighbours_g(i,rose(6)) = -1
536               
537             ENDIF
538          ENDIF
539
540       ENDDO
541
542       IF ( bavard .GT. 1 ) THEN
543          WRITE(numout,*) '  > Total number of points: ',npts_glo
544          WRITE(numout,*) '  > Using default zonal resolution at',ndefault_lon,' points.'
545          WRITE(numout,*) '  > Using default meridional resolution at',ndefault_lat,' points.'
546       ENDIF
547       !
548    ENDIF ! (root_prc)
549
550    CALL scatter(neighbours_g,neighbours)
551    CALL scatter(resolution_g,resolution)
552    CALL scatter(area_g,area)
553    CALL bcast(min_resol)
554    CALL bcast(max_resol)
555    IF ( bavard .EQ. 5 ) THEN
556       WRITE(numout,*) '  > resolution  = ',resolution
557       WRITE(numout,*) '  > rose = ',rose
558       WRITE(numout,*) '  > neighbours  = ',neighbours
559    ENDIF
560    IF ( bavard .GT. 1 ) WRITE(numout,*) 'Leaving grid_stuff'
561   
562  END SUBROUTINE grid_stuff
563  !
564!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
565
566END MODULE grid
Note: See TracBrowser for help on using the repository browser.