--- MITgcm/pkg/fizhi/fizhi_init_veg.F 2004/06/09 20:33:37 1.9 +++ MITgcm/pkg/fizhi/fizhi_init_veg.F 2004/07/23 22:32:28 1.15 @@ -1,8 +1,11 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/fizhi/fizhi_init_veg.F,v 1.9 2004/06/09 20:33:37 molod Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/fizhi/fizhi_init_veg.F,v 1.15 2004/07/23 22:32:28 molod Exp $ C $Name: $ +#include "FIZHI_OPTIONS.h" + subroutine fizhi_init_veg(mythid,vegdata,im,jm,Nsx,Nsy,Nxg,Nyg, - .maxtyp,nchp,lons,lats,surftype,tilefrac,igrd,ityp,chfr,chlt,chlon) + . 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 @@ -51,15 +54,15 @@ implicit none #include "EEPARAMS.h" - integer mythid,im,jm,maxtyp,nchp,Nsx,Nsy,Nxg,Nyg + integer mythid,im,jm,maxtyp,nchp,nchptot,nchpland,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 lats(im,jm,nSx,nSy), lons(im,jm,nSx,nSy) - real chfr(nchp,Nsx,Nsy),chlt(nchp,Nsx,Nsy),chlon(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 nchplocal,nchpland,biglobal,bjglobal + integer biglobal,bjglobal integer*4 im_32, jm_32, Nxg_32, Nyg_32 integer*4 iveg_32(im,jm,maxtyp,Nxg,Nyg) @@ -67,6 +70,14 @@ 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', @@ -80,6 +91,7 @@ call exit(101) endif close(kveg) + IF (myThid.eq.1) THEN imdata = im_32 jmdata = jm_32 @@ -95,13 +107,16 @@ 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 @@ -111,61 +126,92 @@ enddo enddo -c create chip arrays for : + ENDDO + ENDDO + + ENDIF +#endif + +#ifdef ALLOW_MNC + _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 + +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 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) - chlon(nchplnd,bi,bj) = lons(i,j,bi,bj) - chlt (nchplnd,bi,bj) = lats(i,j,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) - chlon(nchplocal,bi,bj) = lons(i,j,bi,bj) - chlt (nchplocal,bi,bj) = lats(i,j,bi,bj) - endif - enddo - enddo - enddo +c nchpland<=nchptot is the actual number of land chips - print *, 'bi ',bi,' bj ',bj - print *, 'Number of Total Tiles: ',nchplocal - print *, 'Number of Land Tiles: ',nchplnd - print * + DO BJ = myByLo(myThid), myByHi(myThid) + DO BI = myBxLo(myThid), myBxHi(myThid) +c land points +c ----------- + nchpland = 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 = nchpland + 1 + igrd (nchpland,bi,bj) = i + (j-1)*im + ityp (nchpland,bi,bj) = surftype(i,j,k,bi,bj) + chfr (nchpland,bi,bj) = tilefrac(i,j,k,bi,bj) + chlon(nchpland,bi,bj) = lons(i,j,bi,bj) + chlt (nchpland,bi,bj) = lats(i,j,bi,bj) + endif + enddo + enddo + enddo + +c ocean points +c ------------ + nchptot = nchpland + + 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 = nchptot + 1 + igrd (nchptot,bi,bj) = i + (j-1)*im + ityp (nchptot,bi,bj) = surftype(i,j,k,bi,bj) + chfr (nchptot,bi,bj) = tilefrac(i,j,k,bi,bj) + chlon(nchptot,bi,bj) = lons(i,j,bi,bj) + chlt (nchptot,bi,bj) = lats(i,j,bi,bj) + endif + enddo + enddo + enddo + + if(bi.eq.1.and.bj.eq.1)then + print *, 'Number of Total Tiles: ',nchptot + print *, 'Number of Land Tiles: ',nchpland + print * + endif + + ENDDO ENDDO - ENDDO - + RETURN END