source: branches/publications/ORCHIDEE_CAN_r3069/src_global/grid.f90 @ 7189

Last change on this file since 7189 was 2945, checked in by sebastiaan.luyssaert, 9 years ago

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