pkg/fizhi/fizhi_init_veg.F

C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_init_veg.F,v 1.5 2004/06/07 20:26:43 molod Exp $
C $Name:  $

      subroutine fizhi_init_veg(mythid,vegdata,im,jm,Nsx,Nsy,Nxg,Nyg,
     . maxtyp,nchp,surftype,tilefrac,igrd,ityp,chfr )
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       - model grid longitude dimension
C jm       - model grid 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
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"

      integer mythid,im,jm,maxtyp,nchp,Nsx,Nsy,Nxg,Nyg
      integer surftype(im,jm,maxtyp,Nsx,Nsy)
      integer igrd(nchp,Nsx,Nsy),ityp(nchp,Nsx,Nsy)
      real tilefrac(im,jm,maxtyp,Nsx,Nsy)
      real chfr(nchp,Nsx,Nsy)
      character*40 vegdata
      integer imdata,jmdata,Nxgdata,Nygdata
      integer nchplocal,nchpland,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

      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
                                                                                 
      DO BJ = myByLo(myThid), myByHi(myThid)
      DO BI = myBxLo(myThid), myBxHi(myThid)

      biglobal=bi+(myXGlobalLo-1)/im
      bjglobal=bj+(myYGlobalLo-1)/jm
                                                                                 
      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

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

c  nchplnd<=nchplocal is the actual number of land chips

c land points
c -----------
      nchplnd = 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
             nchplnd  = nchplnd + 1
       igrd (nchplnd,bi,bj) = i + (j-1)*im
       ityp (nchplnd,bi,bj) = surftype(i,j,k,bi,bj)
       chfr (nchplnd,bi,bj) = tilefrac(i,j,k,bi,bj)
      endif
      enddo
      enddo
      enddo

c ocean points
c ------------
      nchplocal = nchplnd

      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
             nchplocal  = nchplocal + 1
       igrd (nchplocal,bi,bj) = i + (j-1)*im
       ityp (nchplocal,bi,bj) = surftype(i,j,k,bi,bj)
       chfr (nchplocal,bi,bj) = tilefrac(i,j,k,bi,bj)
      endif
      enddo
      enddo
      enddo

      print *, 'bi ',bi,' bj ',bj
      print *, 'Number of Total Tiles: ',nchplocal
      print *, 'Number of Land  Tiles: ',nchplnd
      print *

      ENDDO
      ENDDO

      RETURN
      END
1	molod	1.6	C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_init_veg.F,v 1.5 2004/06/07 20:26:43 molod Exp $
2	molod	1.1	C $Name: $
3
4	molod	1.6	subroutine fizhi_init_veg(mythid,vegdata,im,jm,Nsx,Nsy,Nxg,Nyg,
5			. maxtyp,nchp,surftype,tilefrac,igrd,ityp,chfr )
6	molod	1.1	C***********************************************************************
7	molod	1.2	C Subroutine fizhi_init_veg - routine to read in the land surface types,
8			C interpolate to the models grid, and set up tile space for use by
9			C the land surface model, the albedo calculation and the surface
10			C roughness calculation.
11	molod	1.1	C
12			C INPUT:
13			C
14			C mythid - thread number (processor number)
15	molod	1.2	C vegdata - Character*40 Vegetation Dataset name
16	molod	1.1	C im - model grid longitude dimension
17			C jm - model grid latitude dimension (number of lat. points)
18	molod	1.4	C Nsx - Number of processors in x-direction
19			C Nsy - Number of processors in y-direction
20	molod	1.1	C maxtyp - maximum allowable number of land surface types per grid box
21	molod	1.6	C nchp - integer per-processor number of tiles in tile space
22	molod	1.1	C
23			C OUTPUT:
24			C
25	molod	1.4	C surftype - integer array of land surface types [im,jm,maxtyp,Nsx,Nsy]
26	molod	1.1	C tilefrac - real array of corresponding land surface type fractions
27	molod	1.4	C [im,jm,maxtyp,Nsx,Nsy]
28	molod	1.1	C igrd - integer array in tile space of grid point number for each
29	molod	1.4	C tile [nchp,Nsx,Nsy]
30	molod	1.1	C ityp - integer array in tile space of land surface type for each
31	molod	1.4	C tile [nchp,Nsx,Nsy]
32	molod	1.1	C chfr - real array in tile space of land surface type fraction for
33	molod	1.4	C each tile [nchp,Nsx,Nsy]
34	molod	1.1	C
35			C NOTES:
36			C Vegetation type as follows:
37			C 1: BROADLEAF EVERGREEN TREES
38			C 2: BROADLEAF DECIDUOUS TREES
39			C 3: NEEDLELEAF TREES
40			C 4: GROUND COVER
41			C 5: BROADLEAF SHRUBS
42			C 6: DWARF TREES (TUNDRA)
43			C 7: BARE SOIL
44			C 8: DESERT
45			C 9: GLACIER
46			C 10: DARK DESERT
47			C 100: OCEAN
48			C***********************************************************************
49			implicit none
50	molod	1.5	#include "EEPARAMS.h"
51	molod	1.1
52	molod	1.6	integer mythid,im,jm,maxtyp,nchp,Nsx,Nsy,Nxg,Nyg
53	molod	1.4	integer surftype(im,jm,maxtyp,Nsx,Nsy)
54	molod	1.6	integer igrd(nchp,Nsx,Nsy),ityp(nchp,Nsx,Nsy)
55	molod	1.4	real tilefrac(im,jm,maxtyp,Nsx,Nsy)
56	molod	1.6	real chfr(nchp,Nsx,Nsy)
57	molod	1.1	character*40 vegdata
58	molod	1.6	integer imdata,jmdata,Nxgdata,Nygdata
59			integer nchplocal,nchpland,biglobal,bjglobal
60	molod	1.1
61	molod	1.6	integer*4 im_32, jm_32, Nxg_32, Nyg_32
62			integer*4 iveg_32(im,jm,maxtyp,Nxg,Nyg)
63			real*4 veg_32(im,jm,maxtyp,Nxg,Nyg)
64	molod	1.1
65	molod	1.3	integer i,j,k,bi,bj,ierr1,kveg
66	molod	1.1
67			call mdsfindunit( kveg, myThid )
68	molod	1.2	close(kveg)
69			open(kveg,file=vegdata,form='unformatted',access='sequential',
70			. iostat=ierr1)
71	molod	1.1	if( ierr1.eq.0 ) then
72	molod	1.6	read(kveg)im_32,jm_32,Nxg_32,Nyg_32,IVEG_32,VEG_32
73	molod	1.1	else
74			print *
75			print *, 'Veg Dataset: ',vegdata,' not found!'
76			print *
77	molod	1.2	call exit(101)
78	molod	1.1	endif
79	molod	1.2	close(kveg)
80	molod	1.3	IF (myThid.eq.1) THEN
81	molod	1.2	imdata = im_32
82			jmdata = jm_32
83	molod	1.6	Nxgdata = Nxg_32
84			Nygdata = Nyg_32
85	molod	1.2	if( (imdata.ne.im) .or. (jmdata.ne.jm) .or.
86	molod	1.6	. (Nxgdata.ne.Nxg) .or. (Nygdata.ne.Nyg) ) then
87	molod	1.2	print *
88			print *, 'Veg Data Resolution is Incorrect! '
89			print *,' Model Res: ',im,'x',jm,' Data Res: ',imdata,'x',jmdata
90	molod	1.6	print *,' Model Nxg Nyg: ',Nxg,' ',Nyg,' Data Nxg Nyg: ',Nxgdata,
91			. ' ',Nygdata
92	molod	1.2	print *
93			call exit(102)
94	molod	1.1	ENDIF
95	molod	1.2
96	molod	1.5	DO BJ = myByLo(myThid), myByHi(myThid)
97			DO BI = myBxLo(myThid), myBxHi(myThid)
98
99			biglobal=bi+(myXGlobalLo-1)/im
100			bjglobal=bj+(myYGlobalLo-1)/jm
101	molod	1.2
102			do k = 1,maxtyp
103	molod	1.4	do j = 1,jm
104			do i = 1,im
105	molod	1.5	surftype(i,j,k,bi,bj) = iveg_32(i,j,k,biglobal,bjglobal)
106			tilefrac(i,j,k,bi,bj) = veg_32(i,j,k,biglobal,bjglobal)
107	molod	1.2	enddo
108	molod	1.1	enddo
109			enddo
110
111			c create chip arrays for :
112			c igrd : grid index
113			c ityp : veg. type
114			c chfr : vegetation fraction
115
116	molod	1.6	c nchplnd<=nchplocal is the actual number of land chips
117	molod	1.1
118			c land points
119			c -----------
120			nchplnd = 0
121			do k=1,maxtyp
122			do j=1,jm
123			do i=1,im
124	molod	1.2	if(surftype(i,j,k,bi,bj).lt.100.and.
125			. tilefrac(i,j,k,bi,bj).gt.0.)then
126	molod	1.1	nchplnd = nchplnd + 1
127	molod	1.2	igrd (nchplnd,bi,bj) = i + (j-1)*im
128			ityp (nchplnd,bi,bj) = surftype(i,j,k,bi,bj)
129			chfr (nchplnd,bi,bj) = tilefrac(i,j,k,bi,bj)
130	molod	1.1	endif
131			enddo
132			enddo
133			enddo
134
135			c ocean points
136			c ------------
137	molod	1.6	nchplocal = nchplnd
138	molod	1.1
139			do k=1,maxtyp
140			do j=1,jm
141			do i=1,im
142	molod	1.2	if(surftype(i,j,k,bi,bj).ge.100 .and.
143			. tilefrac(i,j,k,bi,bj).gt.0.)then
144	molod	1.6	nchplocal = nchplocal + 1
145			igrd (nchplocal,bi,bj) = i + (j-1)*im
146			ityp (nchplocal,bi,bj) = surftype(i,j,k,bi,bj)
147			chfr (nchplocal,bi,bj) = tilefrac(i,j,k,bi,bj)
148	molod	1.1	endif
149			enddo
150			enddo
151			enddo
152
153	molod	1.3	print *, 'bi ',bi,' bj ',bj
154	molod	1.6	print *, 'Number of Total Tiles: ',nchplocal
155	molod	1.1	print *, 'Number of Land Tiles: ',nchplnd
156			print *
157	molod	1.2
158			ENDDO
159			ENDDO
160	molod	1.1
161			RETURN
162			END