SIB: balan.F90 Source File

Go to the documentation of this file.
00001 !...end-of-timestep calculation to determine if water and energy
00002 !...balance is maintained
00003 
00004 subroutine balan(sib,nsl_old,wwwliq_old,wwwice_old,cas_q)
00005 
00006 use kinds
00007 use sibtype
00008 use sib_const_module, only: &
00009     nsoil,  &
00010     snofac, &
00011     dtt,    &
00012     dti,    &
00013     tau, &
00014     denh2o, &
00015     denice
00016 
00017 use physical_parameters, only: &
00018     hltm
00019 
00020 use eau_params, only : &
00021     lfus
00022 
00023 
00024 implicit none
00025 
00026 !-----------------------------------------------------------------------
00027 !...WATER BALANCE...
00028 !
00029 !     WATER IN + CHANGE IN STORAGE = WATER OUT
00030 !
00031 !     inputs:
00032 !       precipitation             (sib_prog%cupr,sib_prog%lspr)    (mm/sec)
00033 !
00034 !     output2: 
00035 !       CAS-mixed layer latent heat flux   (sib_diag%fws)          (W/m^2)
00036 !
00037 !       runoff
00038 !          surface runoff                   (sib_diag%roffo)       (mm/hr)
00039 !          subsurface runoff                (sib_diag%roff)        (mm/hr)
00040 !
00041 !       storage
00042 !          interception reservoirs
00043 !             canopy interception           (sib_prog%capac(1))    (kg/m^2)
00044 !             ground interception           (sib_prog%capac(2))    (kg/m^2)
00045 !          soil
00046 !             soil water                    (sib_prog%www_liq)     (kg/m^2)
00047 !             soil ice                      (sib_prog%www_ice)     (kg/m^2)
00048 !          snow
00049 !             snow water                    (sib_prog%www_liq)     (kg/m^2)
00050 !             snow ice                      (sib_prog%www_ice)     (kg/m^2)
00051 !
00052 !           canopy airspace                  (cas_q)           (kg/m^2)
00053 !
00054 !              ** note: indices for soil water and ice arrays
00055 !                       will be 1 to 10 for soil (1 at top,
00056 !                       10 at bottom of soil column), and from
00057 !                       -4 (top) to 0 (bottom) for up to 5 
00058 !                       snow layers.
00059 !
00060 !----------------------------------------------------------------------
00061 !...ENERGY BALANCE
00062 !
00063 !    ENERGY IN + CHANGE IN STORAGE = ENERGY OUT
00064 !
00065 !    inputs:
00066 !      radiation absorbed by canopy    (sib_diag%radt(1))       (W/m^2)
00067 !      radiation absorbed by ground    (sib_diag%radt(2))       (W/m^2)
00068 !
00069 !    outputs:   
00070 !      latent heat 
00071 !       transpiration                       (sib_diag%ect)      (J/m^2)
00072 !       evaporation
00073 !          from canopy interception storage (sib_diag%eci)      (J/m^2)
00074 !          from ground interception storage (sib_diag%egi)      (J/m^2)
00075 !          from surface soil layer          (sib_diag%egs)      (J/m^2)
00076 !          from snow                        (sib_diag%ess)      (J/m^2)  
00077 !      sensible heat
00078 !       canopy                              (sib_diag%hc)       (J/m^2)
00079 !       ground                              (sib_diag%hg)       (J/m^2)
00080 !       snow                                (sib_diag%hs)       (J/m^2)
00081 !
00082 !    storage
00083 !      canopy storage flux                  (sib_diag%chf)      (W/m^2)
00084 !      soil storage flux                    (sib_diag%shf)      (W/m^2)
00085 !
00086 !
00087 !----------------------------------------------------------------------
00088 
00089 type(sib_t), intent(inout) :: sib
00090 
00091 !-----------------------------------------------------------------------
00092 !     input variables
00093 !-----------------------------------------------------------------------
00094 
00095 integer(kind=int_kind),intent(in)      :: nsl_old
00096 real(kind=dbl_kind),intent(in), 
00097     dimension(-nsnow+1:nsoil) :: wwwliq_old    ! (kg/m^2)
00098 real(kind=dbl_kind),intent(in), 
00099     dimension(-nsnow+1:nsoil) :: wwwice_old    ! (kg/m^2)
00100 
00101 real(kind=dbl_kind),intent(in) :: cas_q         ! (kg/m^2)
00102 
00103 !-----------------------------------------------------------------------
00104 !     end input variables
00105 !-----------------------------------------------------------------------
00106 !     local variables
00107 !-----------------------------------------------------------------------
00108 
00109 real(kind=dbl_kind) :: dstor    ! change in storage, canopy and 
00110 !   surface interception (kg/m^2)
00111 
00112 real(kind=dbl_kind)   :: dqsoil   ! change in storage, soil only
00113 !   liquid and ice (kg/m^2)
00114 
00115 real(kind=dbl_kind) :: dqsnow   ! change in storage, snow (kg/m^2)
00116 real(kind=dbl_kind) :: dqvegsnow   ! change in storage, snow (kg/m^2)
00117 real(kind=dbl_kind) :: snownew  ! amount of snow, end-of-timestep (kg/m^2)
00118 real(kind=dbl_kind) :: snowold  ! amount of snow, beg-of-timestep (kg/m^2)
00119 
00120 real(kind=dbl_kind) :: evap     ! evaporation, from CAS to atmosphere.
00121 !   this encompasses CAS storage, 
00122 !   ground interception stores, and
00123 !   evap from top soil layer as well 
00124 !   as snow (kg/m^2)
00125 
00126 real(kind=dbl_kind) :: transp   ! transpiration (kg/m^2)
00127 
00128 real(kind=dbl_kind) :: runoff   ! total runoff: surface + subsurface
00129 !    (kg/m^2)
00130 real(kind=dbl_kind) :: sbeg,send! sum holders
00131 
00132 real(kind=dbl_kind) :: precip   ! precip (kg/m^2)
00133 
00134 real(kind=dbl_kind) :: cas_q_new ! CAS water stored as vapor (kg/m^2)
00135 
00136 real(kind=dbl_kind) :: rhs,lhs  ! right and left hand sides of energy
00137 !  balance equation (W/m^2)
00138 
00139 integer(kind=int_kind) :: i        ! loop index
00140 
00141 
00142 
00143     !-----------------------------------------------------------------------
00144     !     end local variables
00145     !-----------------------------------------------------------------------
00146 
00147     !...NOTE...unless otherwise noted, all units will be kg/m^2 water
00148 
00149 
00150     !...change in canopy and surface interception storage 
00151 ! CSR moved capac_old to diag
00152     dstor = (sib%prog%capac(1) - sib%diag%capac_old(1))   +  & ! canopy interception
00153         (sib%prog%capac(2) - sib%diag%capac_old(2))        ! surface interception
00154 
00155     !...change in soil water
00156     dqsoil = 0.0
00157     do i=1,nsoil
00158         dqsoil = dqsoil + (sib%prog%www_liq(i) - wwwliq_old(i)) + &
00159             (sib%prog%www_ice(i) - wwwice_old(i))
00160     enddo
00161 
00162     !...beginning of timestep snow water (liquid + ice)
00163 ! CSR count full snow column
00164     snowold = 0.0
00165     do i=-nsnow+1,0
00166         snowold = snowold + wwwliq_old(i) + wwwice_old(i)
00167     enddo
00168 
00169     !...end of timestep snow water (liquid + ice)
00170 ! CSR count full snow column
00171     snownew = 0.0
00172     do i=-nsnow+1,0
00173         snownew = snownew + sib%prog%www_liq(i) + sib%prog%www_ice(i)
00174     enddo
00175 
00176 !    dqsnow = snownew - snowold
00177     dqsnow = sib%prog%snow_mass - sib%diag%snow_mass_old
00178     dqvegsnow = sib%prog%snow_veg - sib%diag%snow_veg_old
00179 
00180     !...evaporation, from canopy and ground surface interception 
00181     !...stores (sib%diag%eci, sib%diag%egi), 
00182     !...as well as top soil layer (sib%diag%egs) and snow evaporation
00183 
00184     evap = sib%diag%fws * dtt / hltm           ! total CAS water vapor flux
00185 
00186     !...runoff, surface + subsurface
00187     runoff = (sib%diag%roffo + sib%diag%roff)
00188 
00189     !...precip
00190     precip = (sib%prog%lspr + sib%prog%cupr) * dtt
00191 
00192     !...WATER BALANCE...
00193 ! CSR: works pretty well now
00194     sib%diag%wbal =  precip - (evap + runoff) - &
00195         (dqsoil + dqsnow + dqvegsnow + dstor) - sib%diag%cas_w_storage * dtt / hltm
00196 
00197 
00198     !...setting water balance error at 0.1 kg/m^2 (or 0.1mm)
00199     if(abs(sib%diag%wbal) > 0.5_dbl_kind) then
00200         print*,'WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW'
00201         print'(2(a,g16.6))','WATER IMBALANCE: hour=',tau,'imbalance=' &
00202             ,sib%diag%wbal,' kg/m^2'
00203         print'(a)','TERMS:  input - output - storage (all units kg/m^2)'
00204         print'(3(a,g14.6))','input=',precip,' output=',evap+runoff &
00205             ,' storage=',dqsoil + dqsnow + dstor + dqvegsnow + sib%diag%cas_w_storage * dtt / hltm
00206         print'(2(a,g16.6))','evap(fws)=',evap,' runoff=',runoff
00207         print*,' '
00208         print'(a,4g14.6)','output (eci,egi,egs,ess):',sib%diag%eci/hltm &
00209             ,sib%diag%egi/hltm,sib%diag%egs/hltm &
00210             ,sib%diag%ess*snofac/hltm
00211         print'(a,5f14.6)','storage (interception, soil, snow,vapor, canopy snow):' &
00212             ,dstor,dqsoil,dqsnow,sib%diag%cas_w_storage * dtt / hltm, dqvegsnow
00213         print'(a,2f14.6)','runoff (overland, subsfc):' &
00214             ,sib%diag%roffo  ,sib%diag%roff
00215 
00216         print*,'soil moisture'
00217         sbeg = 0.0
00218         send = 0.0
00219         do i=1,nsoil
00220             print'(i4,2(a,g14.6))',i,&
00221                  ' beg liq=',wwwliq_old(i), ' end liq=',sib%prog%www_liq(i)
00222             sbeg = sbeg + wwwliq_old(i)
00223             send = send + sib%prog%www_liq(i)
00224         enddo
00225         print'(2(a,g14.6))','sum beg=',sbeg,' sum end=',send
00226 
00227         print*,'Dcapac1: ',sib%prog%capac(1) - sib%diag%capac_old(1),sib%prog%capac(1)
00228         print*,'Dcapac2: ',sib%prog%capac(2) - sib%diag%capac_old(2),sib%prog%capac(2)
00229        
00230 
00231         sbeg = 0.0
00232         send = 0.0    
00233         print*,'soil ice'
00234         do i=1,nsoil
00235             print'(i4,2(a,g14.6))',i,' beg ice=',wwwice_old(i),' end ice=' &
00236                 ,sib%prog%www_ice(i)
00237             sbeg = sbeg + wwwice_old(i)
00238             send = send + sib%prog%www_ice(i)
00239         enddo
00240         print'(2(a,g14.6))','sum beg=',sbeg,' sum end=',send
00241 
00242         print*,'snow',nsl_old,sib%prog%nsl,sib%diag%snow_mass_old,sib%prog%snow_mass
00243         sbeg = 0.0
00244         send = 0.0
00245         do i=-nsnow+1,0
00246             print'(i4,2(a,g14.6))',i,' beg ice=',wwwice_old(i),' end ice=' ,sib%prog%www_ice(i)
00247             sbeg = sbeg + wwwice_old(i)
00248             send = send + sib%prog%www_ice(i)
00249         enddo
00250 
00251         print'(2(a,g14.6))','sum beg=',sbeg,' sum end=',send
00252         sbeg = 0.0
00253         send = 0.0
00254         do i=-nsnow+1,0
00255             print'(i4,2(a,g14.6))',i,' beg liq=',wwwliq_old(i),' end liq=' ,sib%prog%www_liq(i)
00256             sbeg = sbeg + wwwliq_old(i)
00257             send = send + sib%prog%www_liq(i)
00258         enddo
00259         print'(2(a,g14.6))','sum beg=',sbeg,' sum end=',send
00260 
00261 !        print*,'WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW'
00262 !        stop
00263     endif
00264 
00265     !----------------- END OF WATER BALANCE CODE ----------------------------
00266 
00267    !----------------- ENERGY BALANCE ---------------------------------------
00268 
00269 ! CSR where has the heat term from melt/freeze energy gone?
00270 
00271 ! CSR canopy air space balance (W/m2) - still has spikes (snow? eci/egi?)
00272 
00273     rhs = sib%diag%cas_w_storage + &  ! change in CAS water storage
00274           sib%diag%cas_e_storage + &    ! change in CAS heat storage
00275           sib%diag%fss      +      &  ! sensible heat flux
00276           sib%diag%fws               ! latent heat flux
00277 
00278     lhs = (sib%diag%hg + sib%diag%hc + sib%diag%hs + &
00279           sib%diag%ect + sib%diag%eci + sib%diag%egi + &
00280           sib%diag%egs  + sib%diag%ess) * dti
00281         
00282     sib%diag%abal = lhs - rhs
00283 
00284 ! CSR canopy balance (W/m2) - still has spikes (snow? eci/egi?)
00285     sib%diag%cbal = sib%diag%radtt(1) - sib%diag%chf - &
00286          ( sib%diag%ect + sib%diag%eci + sib%diag%hc ) * dti
00287 
00288 ! CSR ground balance (W/m2) - way off +/- 50W/m2
00289     sib%diag%gbal = sib%diag%radtt(2) +  sib%diag%radtt(3) & 
00290          - sib%diag%shf &
00291          - (sib%diag%hg + sib%diag%hs + sib%diag%egi + sib%diag%egs + sib%diag%ess ) * dti
00292 
00293 ! CSR total energy balance (W/m2) - way off +/- 50W/m2
00294 
00295     rhs = sib%diag%cas_e_storage + sib%diag%cas_w_storage + &
00296          sib%diag%fss + sib%diag%fws + sib%diag%chf + sib%diag%shf
00297 
00298     lhs = sib%diag%radtt(1) + sib%diag%radtt(2) + sib%diag%radtt(3)
00299 
00300     sib%diag%ebal = lhs - rhs
00301 
00302     if (abs(sib%diag%ebal) > 0.001) then
00303         !   if (ebal /= 0.0) then
00304         print*,'EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE'
00305         print*, 'point',sib%stat%pt_num
00306         print'(4(a,g14.6))','ENERGY IMBALANCE: hour =',tau,'rhs=' &
00307             ,rhs,' lhs=',lhs,' imbalance=',sib%diag%ebal
00308         print*,'imbalance amount=',sib%diag%ebal,' (W/m^2)'
00309         print*,'canopy air space balance=',sib%diag%abal,' (W/m^2)'
00310         print*,'canopy balance=',sib%diag%cbal,' (W/m^2)'
00311         print*,'ground balance=',sib%diag%gbal,' (W/m^2) AREAS:',sib%diag%areas
00312         print'(2(a,g14.6))','sib%diag%ect=',sib%diag%ect*dti,' sib%diag%eci=' &
00313             ,sib%diag%eci*dti
00314         print'(2(a,g14.6))','sib%diag%egi=',sib%diag%egi*dti,' sib%diag%egs=' &
00315             ,sib%diag%egs*dti
00316         print'(2(a,g14.6))','sib%diag%ess=',sib%diag%ess*dti,' sib%diag%hc='  &
00317             ,sib%diag%hc*dti
00318         print'(2(a,g14.6))','sib%diag%hg=',sib%diag%hg*dti,' sib%diag%hs='    &
00319             ,sib%diag%hs*dti
00320         print'(3(a,g14.6))','radt(1)=',sib%diag%radt(1),' radt(2)=',          &
00321             sib%diag%radt(2),' radt(3)=',sib%diag%radt(3)
00322         print'(3(a,g14.6))','radtt(1)=',sib%diag%radtt(1),' radtt(2)=',          &
00323             sib%diag%radtt(2),' radtt(3)=',sib%diag%radtt(3)
00324         print'(2(a,g14.6))','sib%diag%chf=',sib%diag%chf,' sib%diag%shf=',    &
00325             sib%diag%shf
00326         print*,'EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE'
00327         print*,'snow layers: ',-sib%prog%nsl
00328 
00329 !        stop
00330 
00331     endif
00332 
00333 
00334 end subroutine balan