getmet.f90

CLM集合卡曼滤波数据同化算法

   SUBROUTINE getmet (lon ,lat, lun_met) !,loc_yea,loc_day,loc_hou) 

! ======================================================================
!      Source file: getmet.f90
! Original version: Yongjiu Dai, September 15, 1999
!
! Access meteorological data
!
! ======================================================================

   use enkfpar_module ! logical unit number of obs_file is stored in it
   use enkfobs_module ! sat_obs is read into it 

   USE clm2d_module
   IMPLICIT NONE

   integer, INTENT(in) :: &
         lon,       &! atm number of longitudes
         lat,       &! atm number of latitudes
         lun_met     ! logical unit number of atmospheric forcing data 

!* real, INTENT(out), dimension(lon, lat) :: &
!*       pgcmxy,    &! atm bottom level pressure (pa)
!*       psrfxy,    &! surface pressure (pa)
!*       ugcmxy,    &! atm bottom level zonal wind (m/s)
!*       vgcmxy,    &! atm bottom level meridional wind (m/s)
!*       tgcmxy,    &! atm bottom level temperature (kelvin)
!*       qgcmxy,    &! atm btm level specific humidity (kg/kg)
!*       prcxy,     &! atm convective precipitation rate (mm h2o/s)
!*       prlxy,     &! atm large-scale precipitation rate (mm h2o/s)
!*       flwdsxy,   &! atm downward longwave rad onto surface (w/m**2)
!*       solsxy,    &! atm vis direct beam solar rad onto srf (w/m**2)
!*       sollxy,    &! atm nir direct beam solar rad onto srf (w/m**2)
!*       solsdxy,   &! atm vis diffuse solar rad onto srf (w/m**2)
!*       solldxy     ! atm nir diffuse solar rad onto srf(w/m**2)

!* real, INTENT(out), dimension(lon, lat, 3) :: &
!*       zgcmxy      ! atm btm level height above surface (m)

! local
   real  solar,     &! incident solar radiation [W/m2]
         frl,       &! atmospheric infrared (longwave) radiation [W/m2]
         prcp,      &! precipitation [mm/s]
         tm,        &! temperature at reference height [kelvin]
         us,        &! wind component in eastward direction [m/s]
         vs,        &! wind component in northward direction [m/s]
         pres,      &! atmospheric pressure [pa]
         qm,        &! specific humidity at reference height [kg/kg]
       	 f 			 ! relative humidity [pa/pa]
   integer i, j      ! looping index

   integer IYEAR,IMONTH,IDAY,IHOUR
   
   
   !real rnet,twet,esat,esatdT,qsat,qsatdT,us_d 
   !real rrref, ddd, ggg
   character*2 site
   !real time, VPA
   !integer flag, Rflag
   !real calday, alon, alat, cosz, so, a, r, k, pi
   real::E
   integer::loc_yea,loc_day,loc_hou

! ------------------ note ------------------
! the model required the same longitudinal resolution for
! all latitude strip. For the variable longitudinal resolution
! cases, please assign the meteorological values to 
! -999 to the land grids which are not included in the calcultion.
! ----------------------------------------------------------------------
![3] ABRACOS site Reserva Jaru & Fazenda N.S. 
! (obs height forest: zwind=52.5)
! (obs height pasture: zwind=5.4, ztem=5.1)
!      read(lun_met,*)
     do j=1,lat
        do i=1,lon 
     

!      read(lun_met,*) site,iyear,iday,ihour,&
!                      us,tm,f,pres,solar,frl,prcp !solar,rrref,rnet,twet,tm,us,ddd,ggg,prcp
        read(lun_met,*) us,tm,f,pres,solar,frl,prcp 
       print*,'read sucessfully' 
30    format(a2,i4,i3,i4,7(f10.2))
      
! ##############################################################
!# the following code is used to read MODIS-TEMPERATURE 
!#
!#
! #############################################################
        read(lun_obs,*) sat_obs


!      loc_yea=iyear
!	  loc_day=iday
!	  loc_hou=ihour
 
      !if(solar.eq.-6999.00) print*,'solar',iyear,iday,ihour
      !if(rnet.eq.-6999.00) print*,'rnet',iyear,iday,ihour
      !if(rrref.eq.-6999.00) print*,'rrref',iyear,iday,ihour
      !if(tm+273.16.lt.274.) print*,'tm',iyear,iday,ihour
      !if(twet.eq.-6999.00) print*,'twet',iyear,iday,ihour
      !if(tm.eq.-6999.00) print*,'tm',iyear,iday,ihour
      !if(us.eq.-6999.00) print*,'us',iyear,iday,ihour
      if(prcp.eq.-6999.00) print*,'prcp',iyear,iday,ihour
	  if(frl.eq.-6999.00)print*,'frl ',iyear,iday,ihour

 
! for missing data
      if ((prcp.lt.0.0).or.(prcp.eq.(-6999.00)))   prcp  =0.0
      if (solar.lt.0.0)                            solar =0.0
      if (us.lt.0.0)                               us    =1.0
	  !if (frl.eq.-6999.00)                         frl   =200.0
	  
      
 
      !if(rrref.eq.-99.99) rrref=solar*0.122
      !if(rrref.lt.0.0) rrref=0.0
 
      !if(rnet.eq.-99.99) then
        !if(site.eq.'RD' .or. site.eq.'rd') rnet = 0.79*solar - 26.08
        !if(site.eq.'FD' .or. site.eq.'fd') rnet = 0.79*solar - 16.26
      !endif
 
      !pres =1000.0*100.0 
! specific humidity
       pres=pres*10.0
      if (tm.gt.0.0) then

	  E=6.1078*10**(7.5*tm/(237.3+tm))

          qm=(0.622*f*E/100.0)/(pres-f*E/100.0)

      else 

          E=6.1078*10**(9.5*tm/(265.5+tm))

          qm=(0.622*f*E/100.0)/(pres-f*E/100.0)

      end if
	  
	  
	  vs=0.0


      tm = tm+273.16 
      !if(tm.lt.274.) stop 'meterological data error 2'
      !twet  = twet + 273.16
      !call qsadv(twet,pres,esat,esatdT,qsat,qsatdT)
      !qm = qsat - 3.5*2.8704E2*(tm-twet)/2.50036e6

	  !qm=(qm/100.0)*0.622*((61.1*7.5*(tm-273.16))/(237.3+(tm-273.16)))/(pres-(61.1*7.5*(tm-273.16))/(237.3+(tm-273.16)))
 
      !frl = rnet
      pres=pres*100.0
      prcp=prcp/600.0
! ----------------------------------------------------------------------
!      do j = 1, lat
!         do i = 1, lon
            pgcmxy(i,j)  = pres   
            psrfxy(i,j)  = pres 
            ugcmxy(i,j)  = us   
            vgcmxy(i,j)  = vs  
            tgcmxy(i,j)  = tm 
            qgcmxy(i,j)  = qm  
            prcxy(i,j)   = 0. 
            prlxy(i,j)   = prcp

            flwdsxy(i,j) = frl  
            print*,frl
! relation between fraction diffuse and atmospheric transmission 
! recommended by de Jong (Spitters et al. 1986) for hourly radiation values
!
!           alon = -61.916
!           alat = -10.083
!           pi = 3.141592654
!           calday = float(iday) + float(ihour*60)/86400.
!           call clmzen(1, calday, alon, alat, cosz)
!           so = 1370. *(1. + 0.033*cos(2.*pi*float(iday)/365.))*cosz
!           a = solar / so
!           r = 0.847 - 1.61*cosz + 1.04*cosz**2
!           k = (1.47 - r)/1.66
!           if(a <= 0.22)then
!              solsxy(i,j)  = 0.
!              sollxy(i,j)  = 0.
!              solsdxy(i,j) = 0.5*solar
!              solldxy(i,j) = 0.5*solar
!           else if(a <= 0.35)then
!              solsxy(i,j)  = 0.5*solar*6.4*(a-0.22)**2
!              sollxy(i,j)  = 0.5*solar*6.4*(a-0.22)**2
!              solsdxy(i,j) = 0.5*solar*(1.-6.4*(a-0.22)**2)
!              solldxy(i,j) = 0.5*solar*(1.-6.4*(a-0.22)**2)
!           else if(a <= k)then
!              solsxy(i,j)  = 0.5*solar*(-0.47+1.66*a) 
!              sollxy(i,j)  = 0.5*solar*(-0.47+1.66*a)
!              solsdxy(i,j) = 0.5*solar*(1.47-1.66*a)
!              solldxy(i,j) = 0.5*solar*(1.47-1.66*a)
!           else
!              solsxy(i,j)  = 0.5*solar*(1.-r)
!              sollxy(i,j)  = 0.5*solar*(1.-r)
!              solsdxy(i,j) = 0.5*solar*r
!              solldxy(i,j) = 0.5*solar*r
!           endif
!
!           if(a <= 0.07)then
!              solsxy(i,j)  = 0.
!              sollxy(i,j)  = 0.
!              solsdxy(i,j) = 0.5*solar
!              solldxy(i,j) = 0.5*solar
!           else if(a <= 0.35)then
!              solsxy(i,j)  = 0.5*solar*2.3*(a-0.07)**2
!              sollxy(i,j)  = 0.5*solar*2.3*(a-0.07)**2
!              solsdxy(i,j) = 0.5*solar*(1.-2.3*(a-0.07)**2)
!              solldxy(i,j) = 0.5*solar*(1.-2.3*(a-0.07)**2)
!           else if(a < 0.75)then
!              solsxy(i,j)  = 0.5*solar*(-0.33+1.46*a)
!              sollxy(i,j)  = 0.5*solar*(-0.33+1.46*a)
!              solsdxy(i,j) = 0.5*solar*(1.33-1.46*a)
!              solldxy(i,j) = 0.5*solar*(1.33-1.46*a)
!           else
!              solsxy(i,j)  = 0.5*solar*0.77
!              sollxy(i,j)  = 0.5*solar*0.77
!              solsdxy(i,j) = 0.5*solar*0.23
!              solldxy(i,j) = 0.5*solar*0.23
!           endif
!           write(113,100) 2.*solsxy(i,j), 2.*solsdxy(i,j), solar
         
            solsxy(i,j)  = solar*45./100.
            sollxy(i,j)  = solar*45./100.
            solsdxy(i,j) = solar*5./100.
            solldxy(i,j) = solar*5./100.

            if(prcxy(i,j)+prlxy(i,j) .gt. 0.)then
            solsxy(i,j)  = solar* 0./100.
            sollxy(i,j)  = solar* 0./100.
            solsdxy(i,j) = solar*50./100.
            solldxy(i,j) = solar*50./100.
            endif

            zgcmxy(i,j,1) = 3.5
            zgcmxy(i,j,2) = 3.5
            zgcmxy(i,j,3) = 3.5

           pco2xy(i,j) =  pgcmxy(i,j)*355.e-06
           po2xy(i,j)  =  pgcmxy(i,j)*0.209
         enddo 
      enddo 

!100   format(3f12.3)


      END SUBROUTINE getmet