pkg/fizhi/fizhi_init_veg.F

C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_init_veg.F,v 1.17 2004/10/10 06:08:48 edhill Exp $
C $Name:  $

#include "FIZHI_OPTIONS.h"

      subroutine fizhi_init_veg(mythid,vegdata,im,jm,Nsx,Nsy,Nxg,Nyg,
     . maxtyp,nchp,nchptot,nchpland,lons,lats,surftype,tilefrac,
     . igrd,ityp,chfr,chlt,chlon)
C***********************************************************************
C Subroutine fizhi_init_veg - routine to read in the land surface types,
C      interpolate to the models grid, and set up tile space for use by 
C      the land surface model, the albedo calculation and the surface 
C      roughness calculation.
C
C INPUT:
C 
C mythid   - thread number (processor number)
C vegdata  - Character*40 Vegetation Dataset name
C im       - longitude dimension
C jm       - latitude dimension (number of lat. points)
C Nsx      - Number of processors in x-direction
C Nsy      - Number of processors in y-direction
C maxtyp   - maximum allowable number of land surface types per grid box
C nchp     - integer per-processor number of tiles in tile space
C lons     - longitude in degrees [im,jm,nSx,nSy]
C lats     - latitude in degrees [im,jm,nSx,nSy]
C
C OUTPUT:
C
C surftype - integer array of land surface types [im,jm,maxtyp,Nsx,Nsy]
C tilefrac - real array of corresponding land surface type fractions 
C            [im,jm,maxtyp,Nsx,Nsy]
C igrd     - integer array in tile space of grid point number for each 
C            tile [nchp,Nsx,Nsy]
C ityp     - integer array in tile space of land surface type for each 
C            tile [nchp,Nsx,Nsy]
C chfr     - real array in tile space of land surface type fraction for 
C            each tile [nchp,Nsx,Nsy]
C
C NOTES:
C       Vegetation type as follows:
C                  1:  BROADLEAF EVERGREEN TREES
C                  2:  BROADLEAF DECIDUOUS TREES
C                  3:  NEEDLELEAF TREES
C                  4:  GROUND COVER
C                  5:  BROADLEAF SHRUBS
C                  6:  DWARF TREES (TUNDRA)
C                  7:  BARE SOIL
C                  8:  DESERT    
C                  9:  GLACIER
C                 10:  DARK DESERT
C                100:  OCEAN
C***********************************************************************
      implicit none
#include "EEPARAMS.h"
#include "PARAMS.h"

      integer mythid,im,jm,maxtyp,nchp,Nsx,Nsy,Nxg,Nyg
      integer nchptot(Nsx,Nsy), nchpland(Nsx,Nsy)
      integer surftype(im,jm,maxtyp,Nsx,Nsy)
      integer igrd(nchp,Nsx,Nsy),ityp(nchp,Nsx,Nsy)
      _RL tilefrac(im,jm,maxtyp,Nsx,Nsy)
      _RL lats(im,jm,nSx,nSy), lons(im,jm,nSx,nSy)
      _RL chfr(nchp,Nsx,Nsy),chlt(nchp,Nsx,Nsy),chlon(nchp,Nsx,Nsy)
      character*40 vegdata
      integer imdata,jmdata,Nxgdata,Nygdata
      integer biglobal,bjglobal

      integer*4 im_32, jm_32, Nxg_32, Nyg_32
      integer*4 iveg_32(im,jm,maxtyp,Nxg,Nyg)
         real*4  veg_32(im,jm,maxtyp,Nxg,Nyg)

      integer i,j,k,bi,bj,ierr1,kveg

#ifdef ALLOW_MNC
      character*(MAX_LEN_FNAM) fizhi_veg_bn
#endif

C     Allow for MDSIO format if someday needed
#ifdef ALLOW_MDSIO
      IF ( .FALSE. ) THEN

      call mdsfindunit( kveg, myThid )
      close(kveg)
      open(kveg,file=vegdata,form='unformatted',access='sequential', 
     .                      iostat=ierr1)
      if( ierr1.eq.0 ) then
          read(kveg)im_32,jm_32,Nxg_32,Nyg_32,IVEG_32,VEG_32
      else
       print *
       print *, 'Veg Dataset: ',vegdata,' not found!'
       print *
       call exit(101)
      endif
      close(kveg)

      IF (myThid.eq.1) THEN
      imdata = im_32
      jmdata = jm_32
      Nxgdata = Nxg_32
      Nygdata = Nyg_32
      if( (imdata.ne.im) .or. (jmdata.ne.jm) .or. 
     .                     (Nxgdata.ne.Nxg) .or. (Nygdata.ne.Nyg) ) then
       print *
       print *, 'Veg Data Resolution is Incorrect! '
       print *,' Model Res: ',im,'x',jm,' Data Res: ',imdata,'x',jmdata
       print *,' Model Nxg Nyg: ',Nxg,' ',Nyg,' Data Nxg Nyg: ',Nxgdata,
     .                    ' ',Nygdata
       print *
       call exit(102)
      ENDIF
      ENDIF

      DO BJ = myByLo(myThid), myByHi(myThid)
      DO BI = myBxLo(myThid), myBxHi(myThid)

      biglobal=bi+(myXGlobalLo-1)/im
      bjglobal=bj+(myYGlobalLo-1)/jm
#if defined( _BYTESWAPIO ) && defined( ALLOW_MDSIO )
      call MDS_BYTESWAPR4(im*jm*maxtyp,veg_32(1,1,1,biglobal,bjglobal))
#endif
      do k = 1,maxtyp
      do j = 1,jm
      do i = 1,im
       surftype(i,j,k,bi,bj) = iveg_32(i,j,k,biglobal,bjglobal) 
       tilefrac(i,j,k,bi,bj) = veg_32(i,j,k,biglobal,bjglobal)
      enddo
      enddo
      enddo

      ENDDO
      ENDDO

      ENDIF
#endif

#ifdef ALLOW_MNC
      IF (useMNC) THEN
        _BEGIN_MASTER( myThid )
      
        do i = 1,MAX_LEN_FNAM
          fizhi_veg_bn(i:i) = ' '
        enddo
      
C       The following base name should be handled by some sort of input
C       name parameter in FIZHI_READPARMS() plus a possible size.
      
C       Set the base name     1234567890
        fizhi_veg_bn(1:10) = 'fizhi_veg '
      
        CALL MNC_CW_I_R('I', fizhi_veg_bn, 0,0, 
     &       'surftype',surftype,myThid)
        CALL MNC_CW_RL_R('R', fizhi_veg_bn, 0,0,
     &       'tilefrac', tilefrac, myThid)
      
        _END_MASTER( myThid )
      ENDIF
#endif

c     create chip arrays for :
c      igrd :  grid index
c      ityp :  veg. type
c      chfr :  vegetation fraction
c      chlon:  chip longitude
c      chlt :  chip latitude

c     nchpland<=nchptot is the actual number of land chips

      DO BJ = myByLo(myThid), myByHi(myThid)
        DO BI = myBxLo(myThid), myBxHi(myThid)

c         land points
c         -----------
          nchpland(bi,bj) = 0
          do k=1,maxtyp
            do j=1,jm
              do i=1,im
                if(surftype(i,j,k,bi,bj).lt.100 .and.
     .               tilefrac(i,j,k,bi,bj).gt.0.) then
                  nchpland(bi,bj)  = nchpland(bi,bj) + 1
                  igrd (nchpland(bi,bj),bi,bj) = i + (j-1)*im
                  ityp (nchpland(bi,bj),bi,bj) = surftype(i,j,k,bi,bj)
                  chfr (nchpland(bi,bj),bi,bj) = tilefrac(i,j,k,bi,bj)
                  chlon(nchpland(bi,bj),bi,bj) = lons(i,j,bi,bj)
                  chlt (nchpland(bi,bj),bi,bj) = lats(i,j,bi,bj)
                endif
              enddo
            enddo
          enddo
          
c         ocean points
c         ------------
          nchptot(bi,bj) = nchpland(bi,bj)
          
          do k=1,maxtyp
            do j=1,jm
              do i=1,im
                if(surftype(i,j,k,bi,bj).ge.100 .and. 
     .               tilefrac(i,j,k,bi,bj).gt.0.) then
                  nchptot(bi,bj)  = nchptot(bi,bj) + 1
                  igrd (nchptot(bi,bj),bi,bj) = i + (j-1)*im
                  ityp (nchptot(bi,bj),bi,bj) = surftype(i,j,k,bi,bj)
                  chfr (nchptot(bi,bj),bi,bj) = tilefrac(i,j,k,bi,bj)
                  chlon(nchptot(bi,bj),bi,bj) = lons(i,j,bi,bj)
                  chlt (nchptot(bi,bj),bi,bj) = lats(i,j,bi,bj)
                endif
              enddo
            enddo
          enddo
          
          print *,'Number of Total Tiles for bi=',bi,': ',nchptot(bi,bj)
          print *,'Number of Land  Tiles for bi=',bi,': ',nchpland(bi,bj)

        ENDDO
      ENDDO
      
      RETURN
      END
1	molod	1.18	C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_init_veg.F,v 1.17 2004/10/10 06:08:48 edhill Exp $
2	molod	1.1	C $Name: $
3
4	edhill	1.14	#include "FIZHI_OPTIONS.h"
5
6	molod	1.6	subroutine fizhi_init_veg(mythid,vegdata,im,jm,Nsx,Nsy,Nxg,Nyg,
7	molod	1.15	. maxtyp,nchp,nchptot,nchpland,lons,lats,surftype,tilefrac,
8			. igrd,ityp,chfr,chlt,chlon)
9	molod	1.1	C***********************************************************************
10	molod	1.2	C Subroutine fizhi_init_veg - routine to read in the land surface types,
11			C interpolate to the models grid, and set up tile space for use by
12			C the land surface model, the albedo calculation and the surface
13			C roughness calculation.
14	molod	1.1	C
15			C INPUT:
16			C
17			C mythid - thread number (processor number)
18	molod	1.2	C vegdata - Character*40 Vegetation Dataset name
19	molod	1.7	C im - longitude dimension
20			C jm - latitude dimension (number of lat. points)
21	molod	1.4	C Nsx - Number of processors in x-direction
22			C Nsy - Number of processors in y-direction
23	molod	1.1	C maxtyp - maximum allowable number of land surface types per grid box
24	molod	1.6	C nchp - integer per-processor number of tiles in tile space
25	molod	1.7	C lons - longitude in degrees [im,jm,nSx,nSy]
26			C lats - latitude in degrees [im,jm,nSx,nSy]
27	molod	1.1	C
28			C OUTPUT:
29			C
30	molod	1.4	C surftype - integer array of land surface types [im,jm,maxtyp,Nsx,Nsy]
31	molod	1.1	C tilefrac - real array of corresponding land surface type fractions
32	molod	1.4	C [im,jm,maxtyp,Nsx,Nsy]
33	molod	1.1	C igrd - integer array in tile space of grid point number for each
34	molod	1.4	C tile [nchp,Nsx,Nsy]
35	molod	1.1	C ityp - integer array in tile space of land surface type for each
36	molod	1.4	C tile [nchp,Nsx,Nsy]
37	molod	1.1	C chfr - real array in tile space of land surface type fraction for
38	molod	1.4	C each tile [nchp,Nsx,Nsy]
39	molod	1.1	C
40			C NOTES:
41			C Vegetation type as follows:
42			C 1: BROADLEAF EVERGREEN TREES
43			C 2: BROADLEAF DECIDUOUS TREES
44			C 3: NEEDLELEAF TREES
45			C 4: GROUND COVER
46			C 5: BROADLEAF SHRUBS
47			C 6: DWARF TREES (TUNDRA)
48			C 7: BARE SOIL
49			C 8: DESERT
50			C 9: GLACIER
51			C 10: DARK DESERT
52			C 100: OCEAN
53			C***********************************************************************
54			implicit none
55	molod	1.5	#include "EEPARAMS.h"
56	edhill	1.17	#include "PARAMS.h"
57	molod	1.1
58	molod	1.18	integer mythid,im,jm,maxtyp,nchp,Nsx,Nsy,Nxg,Nyg
59			integer nchptot(Nsx,Nsy), nchpland(Nsx,Nsy)
60	molod	1.4	integer surftype(im,jm,maxtyp,Nsx,Nsy)
61	molod	1.6	integer igrd(nchp,Nsx,Nsy),ityp(nchp,Nsx,Nsy)
62	molod	1.10	_RL tilefrac(im,jm,maxtyp,Nsx,Nsy)
63			_RL lats(im,jm,nSx,nSy), lons(im,jm,nSx,nSy)
64			_RL chfr(nchp,Nsx,Nsy),chlt(nchp,Nsx,Nsy),chlon(nchp,Nsx,Nsy)
65	molod	1.1	character*40 vegdata
66	molod	1.6	integer imdata,jmdata,Nxgdata,Nygdata
67	molod	1.15	integer biglobal,bjglobal
68	molod	1.1
69	molod	1.6	integer*4 im_32, jm_32, Nxg_32, Nyg_32
70			integer*4 iveg_32(im,jm,maxtyp,Nxg,Nyg)
71			real*4 veg_32(im,jm,maxtyp,Nxg,Nyg)
72	molod	1.1
73	molod	1.3	integer i,j,k,bi,bj,ierr1,kveg
74	molod	1.1
75	edhill	1.14	#ifdef ALLOW_MNC
76			character*(MAX_LEN_FNAM) fizhi_veg_bn
77			#endif
78
79			C Allow for MDSIO format if someday needed
80			#ifdef ALLOW_MDSIO
81			IF ( .FALSE. ) THEN
82
83	molod	1.1	call mdsfindunit( kveg, myThid )
84	molod	1.2	close(kveg)
85			open(kveg,file=vegdata,form='unformatted',access='sequential',
86			. iostat=ierr1)
87	molod	1.1	if( ierr1.eq.0 ) then
88	molod	1.6	read(kveg)im_32,jm_32,Nxg_32,Nyg_32,IVEG_32,VEG_32
89	molod	1.1	else
90			print *
91			print *, 'Veg Dataset: ',vegdata,' not found!'
92			print *
93	molod	1.2	call exit(101)
94	molod	1.1	endif
95	molod	1.2	close(kveg)
96	edhill	1.14
97	molod	1.3	IF (myThid.eq.1) THEN
98	molod	1.2	imdata = im_32
99			jmdata = jm_32
100	molod	1.6	Nxgdata = Nxg_32
101			Nygdata = Nyg_32
102	molod	1.2	if( (imdata.ne.im) .or. (jmdata.ne.jm) .or.
103	molod	1.6	. (Nxgdata.ne.Nxg) .or. (Nygdata.ne.Nyg) ) then
104	molod	1.2	print *
105			print *, 'Veg Data Resolution is Incorrect! '
106			print *,' Model Res: ',im,'x',jm,' Data Res: ',imdata,'x',jmdata
107	molod	1.6	print *,' Model Nxg Nyg: ',Nxg,' ',Nyg,' Data Nxg Nyg: ',Nxgdata,
108			. ' ',Nygdata
109	molod	1.2	print *
110			call exit(102)
111	molod	1.1	ENDIF
112	molod	1.12	ENDIF
113
114	molod	1.5	DO BJ = myByLo(myThid), myByHi(myThid)
115			DO BI = myBxLo(myThid), myBxHi(myThid)
116
117			biglobal=bi+(myXGlobalLo-1)/im
118			bjglobal=bj+(myYGlobalLo-1)/jm
119	edhill	1.14	#if defined( _BYTESWAPIO ) && defined( ALLOW_MDSIO )
120	molod	1.13	call MDS_BYTESWAPR4(imjmmaxtyp,veg_32(1,1,1,biglobal,bjglobal))
121			#endif
122	molod	1.2	do k = 1,maxtyp
123	molod	1.4	do j = 1,jm
124			do i = 1,im
125	molod	1.5	surftype(i,j,k,bi,bj) = iveg_32(i,j,k,biglobal,bjglobal)
126			tilefrac(i,j,k,bi,bj) = veg_32(i,j,k,biglobal,bjglobal)
127	molod	1.2	enddo
128	molod	1.1	enddo
129			enddo
130
131	edhill	1.14	ENDDO
132			ENDDO
133
134			ENDIF
135			#endif
136
137			#ifdef ALLOW_MNC
138	edhill	1.16	IF (useMNC) THEN
139			_BEGIN_MASTER( myThid )
140	edhill	1.14
141	edhill	1.16	do i = 1,MAX_LEN_FNAM
142			fizhi_veg_bn(i:i) = ' '
143			enddo
144	edhill	1.14
145	edhill	1.16	C The following base name should be handled by some sort of input
146			C name parameter in FIZHI_READPARMS() plus a possible size.
147	edhill	1.14
148	edhill	1.16	C Set the base name 1234567890
149			fizhi_veg_bn(1:10) = 'fizhi_veg '
150	edhill	1.14
151	edhill	1.16	CALL MNC_CW_I_R('I', fizhi_veg_bn, 0,0,
152			& 'surftype',surftype,myThid)
153			CALL MNC_CW_RL_R('R', fizhi_veg_bn, 0,0,
154			& 'tilefrac', tilefrac, myThid)
155	edhill	1.14
156	edhill	1.16	_END_MASTER( myThid )
157			ENDIF
158	edhill	1.14	#endif
159
160			c create chip arrays for :
161	molod	1.1	c igrd : grid index
162			c ityp : veg. type
163			c chfr : vegetation fraction
164	molod	1.7	c chlon: chip longitude
165			c chlt : chip latitude
166	molod	1.1
167	molod	1.15	c nchpland<=nchptot is the actual number of land chips
168	molod	1.1
169	edhill	1.14	DO BJ = myByLo(myThid), myByHi(myThid)
170			DO BI = myBxLo(myThid), myBxHi(myThid)
171	molod	1.1
172	edhill	1.14	c land points
173			c -----------
174	molod	1.18	nchpland(bi,bj) = 0
175	edhill	1.14	do k=1,maxtyp
176			do j=1,jm
177			do i=1,im
178			if(surftype(i,j,k,bi,bj).lt.100 .and.
179			. tilefrac(i,j,k,bi,bj).gt.0.) then
180	molod	1.18	nchpland(bi,bj) = nchpland(bi,bj) + 1
181			igrd (nchpland(bi,bj),bi,bj) = i + (j-1)*im
182			ityp (nchpland(bi,bj),bi,bj) = surftype(i,j,k,bi,bj)
183			chfr (nchpland(bi,bj),bi,bj) = tilefrac(i,j,k,bi,bj)
184			chlon(nchpland(bi,bj),bi,bj) = lons(i,j,bi,bj)
185			chlt (nchpland(bi,bj),bi,bj) = lats(i,j,bi,bj)
186	edhill	1.14	endif
187			enddo
188			enddo
189			enddo
190
191			c ocean points
192			c ------------
193	molod	1.18	nchptot(bi,bj) = nchpland(bi,bj)
194	edhill	1.14
195			do k=1,maxtyp
196			do j=1,jm
197			do i=1,im
198			if(surftype(i,j,k,bi,bj).ge.100 .and.
199			. tilefrac(i,j,k,bi,bj).gt.0.) then
200	molod	1.18	nchptot(bi,bj) = nchptot(bi,bj) + 1
201			igrd (nchptot(bi,bj),bi,bj) = i + (j-1)*im
202			ityp (nchptot(bi,bj),bi,bj) = surftype(i,j,k,bi,bj)
203			chfr (nchptot(bi,bj),bi,bj) = tilefrac(i,j,k,bi,bj)
204			chlon(nchptot(bi,bj),bi,bj) = lons(i,j,bi,bj)
205			chlt (nchptot(bi,bj),bi,bj) = lats(i,j,bi,bj)
206	edhill	1.14	endif
207			enddo
208			enddo
209			enddo
210
211	molod	1.18	print *,'Number of Total Tiles for bi=',bi,': ',nchptot(bi,bj)
212			print *,'Number of Land Tiles for bi=',bi,': ',nchpland(bi,bj)
213
214	edhill	1.14	ENDDO
215	molod	1.2	ENDDO
216	edhill	1.14
217	molod	1.1	RETURN
218			END