diff --git a/interface/framework/obs_op_pdaf.F90 b/interface/framework/obs_op_pdaf.F90 index 1ce5db81..5966b9fa 100644 --- a/interface/framework/obs_op_pdaf.F90 +++ b/interface/framework/obs_op_pdaf.F90 @@ -71,6 +71,8 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p) #if defined CLMSA USE enkf_clm_mod, & ONLY : clm_varsize, clm_paramarr, clmupdate_swc, clmupdate_T, clmcrns_bd + USE enkf_clm_mod, & + ONLY : clmupdate_lai, clm_begp, clm_endp, clm_patch2gc, clm_patchwt #ifdef CLMFIVE USE clm_instMod, & ONLY : soilstate_inst @@ -128,6 +130,36 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p) ! If no special observation operator is compiled, use point observations lpointobs = .true. +! LAI assimilation assuming statevec contains leafc, slatop, dsladlai for each patch +! two component op to get to LAI obs +! 1 calculate tlai from statevec per patch +! state_p contains leafc of varsize then slatop of varsize and then dsadlai of varsize +! 2 average patches according to weight per gridcell +if (clmupdate_lai==2) then + lpointobs = .false. ! disable general obs_op + ! Note: This fails if there ever is an observation pointing to a gridcell with 0 patches + ! this should not happen, but would cause div by zero. + + do i = 1, dim_obs_p ! loop over all observations + avesm = 0.0 ! use avesm as grid cell average collecter + do j = 1, clm_endp-clm_begp ! loop over all patches + write(*,*) 'DEBUG LAI : obs_index_p(i), clm_patch2gc(j)', obs_index_p(i), clm_patch2gc(j) + if (obs_index_p(i)==clm_patch2gc(j)) then + write(*,*) 'DEBUG LAI : state dsladlai', state_p(j+2*clm_varsize) + if (state_p(j+2*clm_varsize)>0.0) then + avesm = avesm + clm_patchwt(j) * ((state_p(j+1*clm_varsize)*(exp(state_p(j)*state_p(j+2*clm_varsize)) - 1.0))/state_p(j+2*clm_varsize)) ! formula for tlai from leafc + else + avesm = avesm + clm_patchwt(j) *(state_p(j+1*clm_varsize)*state_p(j)) ! 2nd formula + endif ! dsladlai decider + endif ! this patch is in gridcell of observation + write(*,*) 'DEBUG LAI : average ', avesm + enddo ! all patches for obs checked + m_state_p(i) = avesm ! == lai of gridcell where the observation is + enddo ! all observations +endif ! clmupdate_lai == 2 : m_state_p now contains lai for each gridcell with an observation. + + + #if defined CLMSA if (clmupdate_T==1) then diff --git a/interface/model/common/enkf.h b/interface/model/common/enkf.h index f943ad90..1dcd2498 100755 --- a/interface/model/common/enkf.h +++ b/interface/model/common/enkf.h @@ -88,6 +88,8 @@ GLOBAL int nx_local,ny_local,nz_local; GLOBAL int clmupdate_swc; GLOBAL int clmupdate_T; GLOBAL int clmupdate_texture; +GLOBAL int clmupdate_lai; +GLOBAL int clmupdate_lai_params; GLOBAL int clmprint_swc; GLOBAL int clmprint_et; GLOBAL int clmstatevec_allcol; diff --git a/interface/model/common/read_enkfpar.c b/interface/model/common/read_enkfpar.c index b2bfec28..2bf563d2 100755 --- a/interface/model/common/read_enkfpar.c +++ b/interface/model/common/read_enkfpar.c @@ -78,6 +78,8 @@ void read_enkfpar(char *parname) clmupdate_swc = iniparser_getint(pardict,"CLM:update_swc",1); clmupdate_T = iniparser_getint(pardict,"CLM:update_T",0); clmupdate_texture = iniparser_getint(pardict,"CLM:update_texture",0); + clmupdate_lai = iniparser_getint(pardict,"CLM:update_lai",0); + clmupdate_lai_params = iniparser_getint(pardict,"CLM:update_lai_params",0); clmprint_swc = iniparser_getint(pardict,"CLM:print_swc",0); clmprint_et = iniparser_getint(pardict,"CLM:print_et",0); clmstatevec_allcol = iniparser_getint(pardict,"CLM:statevec_allcol",0); diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 82fa382a..a4821e32 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -51,10 +51,14 @@ module enkf_clm_mod ! LST in LST assimilation (clmupdate_T) real(r8),allocatable :: clm_paramarr(:) !hcp CLM parameter vector (f.e. LAI) integer, allocatable :: state_clm2pdaf_p(:,:) !Index of column in hydraulic active state vector (nlevsoi,endc-begc+1) + real(r8),allocatable :: clm_patch2gc(:) ! index array to simplify patch to gridcell averaging + real(r8),allocatable :: clm_patchwt(:) ! weight array to simplify patch to gridcell averaging integer(c_int),bind(C,name="clmupdate_swc") :: clmupdate_swc integer(c_int),bind(C,name="clmupdate_T") :: clmupdate_T ! by hcp integer(c_int),bind(C,name="clmupdate_texture") :: clmupdate_texture integer(c_int),bind(C,name="clmprint_swc") :: clmprint_swc + integer(c_int),bind(C,name="clmupdate_lai") :: clmupdate_lai + integer(c_int),bind(C,name="clmupdate_lai_params") :: clmupdate_lai_params #endif integer(c_int),bind(C,name="clmprint_et") :: clmprint_et integer(c_int),bind(C,name="clmstatevec_allcol") :: clmstatevec_allcol @@ -85,6 +89,10 @@ module enkf_clm_mod ! (currently not used for eclm) logical :: newgridcell !only eclm + ! Arrays for LAI Assimilation + real(r8),allocatable :: tlai(:) + real(r8),allocatable :: old_frac(:) + contains #if defined CLMSA @@ -137,6 +145,50 @@ subroutine define_clm_statevec(mype) clm_statevecsize = clm_statevecsize + 3*((endg-begg+1)*nlevsoi) end if + ! LAI assimilation + if(clmupdate_lai==1) then + ! Allocate array to store total leaf area index per patch + if (allocated(tlai)) deallocate(tlai) + allocate(tlai(clm_endp-clm_begp+1)) + ! Allocate array to store fraction of LAI the patch represents. + if (allocated(old_frac)) deallocate(old_frac) + allocate(old_frac(clm_endp-clm_begp+1)) + + ! #patches values per grid-cell ! clm_varsize not accurate due to variable #patches/grid cell + ! set sizes + clm_varsize = 1 ! each grid cell gets 1 value of LAI + !(clm_endp-clm_begp+1) ! Currently no combination of SWC and LAI DA + + clm_statevecsize = (clm_endg - clm_begg + 1) ! statevector is the size of the num of gridcells + !(clm_endp-clm_begp+1) ! So like this if lai is set it takes priority + + if (clmupdate_lai_params==1) then ! add parameters to the state vector. + ! clm_varsize = (clm_endg - clm_begg+1) ! Currently no combination of SWC and LAI DA + ! clm_statevecsize = (clm_endg - clm_begg+1)*1 ! 1 parameters + clm_varsize = (clm_endp-clm_begp+1)+1 ! Currently no combination of SWC and LAI DA + clm_statevecsize = (clm_endp-clm_begp+1)+1 ! 1 parameters + ! So like this if lai is set it takes priority + end if + + ! if (clmupdate_lai_params==3) then + ! clm_varsize = (clm_endp-clm_begp+1) + ! clm_statevecsize = 3*clm_varsize + ! end if + + end if + + ! for obs_op version allocate additional helper index array + if(clmupdate_lai==2) then + clm_varsize = (clm_endp-clm_begp+1) ! equals number of patches + clm_statevecsize = 3*clm_varsize ! 3 var/params per patch + + if (allocated(clm_patch2gc)) deallocate(clm_patch2gc) + allocate(clm_patch2gc(clm_varsize)) + if (allocated(clm_patchwt)) deallocate(clm_patchwt) + allocate(clm_patchwt(clm_varsize)) + + endif + !hcp LST DA if(clmupdate_T==1) then error stop "Not implemented: clmupdate_T.eq.1" @@ -149,7 +201,7 @@ subroutine define_clm_statevec(mype) #endif IF (allocated(clm_statevec)) deallocate(clm_statevec) - if ((clmupdate_swc/=0) .or. (clmupdate_T/=0) .or. (clmupdate_texture/=0)) then + if ((clmupdate_swc/=0) .or. (clmupdate_T/=0) .or. (clmupdate_texture/=0) .or. (clmupdate_lai/=0)) then !hcp added condition allocate(clm_statevec(clm_statevecsize)) end if @@ -381,11 +433,16 @@ end subroutine cleanup_clm_statevec subroutine set_clm_statevec(tstartcycle, mype) use clm_instMod, only : soilstate_inst, waterstate_inst + use clm_instMod, only : bgc_vegetation_inst use clm_varpar , only : nlevsoi ! use clm_varcon, only: nameg, namec ! use GetGlobalValuesMod, only: GetGlobalWrite use ColumnType , only : col + use PatchType , only : patch + use pftconMod , only : pftcon use shr_kind_mod, only: r8 => shr_kind_r8 + use PhotosynthesisMod, only : params_inst + implicit none integer,intent(in) :: tstartcycle integer,intent(in) :: mype @@ -393,6 +450,7 @@ subroutine set_clm_statevec(tstartcycle, mype) real(r8), pointer :: psand(:,:) real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) + real(r8), pointer :: leafc(:) integer :: i,j,jj,g,c,cc,offset integer :: n_c character (len = 34) :: fn !TSMP-PDAF: function name for state vector output @@ -405,6 +463,8 @@ subroutine set_clm_statevec(tstartcycle, mype) psand => soilstate_inst%cellsand_col pclay => soilstate_inst%cellclay_col porgm => soilstate_inst%cellorg_col + leafc => bgc_vegetation_inst%cnveg_carbonstate_inst%leafc_patch + #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble < 0) THEN @@ -429,6 +489,71 @@ subroutine set_clm_statevec(tstartcycle, mype) error stop "Not implemented clmupdate_swc.eq.2" endif + + ! LAI assimilation + ! Case 1: LAI all patches + ! Case 1: Transformations in Set and Update functions + if(clmupdate_lai==1) then + clm_statevec(:) = 0._r8 ! reset statevector to 0 because we add to it for the average + do i=clm_begp,clm_endp ! iterate through patches + ! update the leaf area index based on leafC and SLA + ! Eq 3 from Thornton and Zimmerman, 2007, J Clim, 20, 3902-3923. + ! (slatop(ivt(p))*(exp(leafc(p)*dsladlai(ivt(p))) - 1._r8))/dsladlai(ivt(p)) + if (pftcon%dsladlai(patch%itype(i)) > 0._r8) then + tlai(i) = (pftcon%slatop(patch%itype(i))*(exp(leafc(i)*pftcon%dsladlai(patch%itype(i))) - 1._r8)) & + /pftcon%dsladlai(patch%itype(i)) + else + tlai(i) = pftcon%slatop(patch%itype(i)) * leafc(i) + end if + tlai(i) = max(0._r8, tlai(i)) ! don't allow negative LAI + ! Add to grid cell average of patch with weighted tlai + clm_statevec(patch%gridcell(i)) = clm_statevec(patch%gridcell(i)) + patch%wtgcell(i)*tlai(i) + end do + + ! Require second loop to calculate fraction from collected grid cell average + ! second loop once grid cell average is known calculate fraction + do i=clm_begp,clm_endp + if (clm_statevec(patch%gridcell(i)) > 0._r8) then ! divide by zero protection + old_frac(i) = patch%wtgcell(i)*tlai(i) / clm_statevec(patch%gridcell(i)) + else + old_frac(i) = 0._r8 + end if + end do + endif + + ! Case 2: statevec with leafc, slatop, dsladlai for transform in obs_op + if (clmupdate_lai==2) then + cc = 1 + do i=clm_begp,clm_endp + clm_statevec(cc) = leafc(i) + clm_statevec(cc + 1*clm_varsize) = pftcon%slatop(patch%itype(i)) + clm_statevec(cc + 2*clm_varsize) = pftcon%dsladlai(patch%itype(i)) + clm_patch2gc(cc) = patch%gridcell(i) + clm_patchwt(cc) = patch%wtgcell(i) + cc = cc + 1 + enddo + write(*,*) 'DEBUG LAI : statevec ', clm_statevec(:) + write(*,*) 'DEBUG LAI : patches ', clm_patch2gc(:), clm_patchwt(:) + endif + + if (clmupdate_lai_params==1) then + cc = 1 + do i=clm_begp,clm_endp + clm_statevec(cc+1*clm_varsize+offset) = pftcon%slatop(patch%itype(i)) + cc = cc + 1 + end do + endif + + ! if (clmupdate_lai_params==3) then + ! cc = 1 + ! do i=clm_begp,clm_endp + ! clm_statevec(cc+1*clm_varsize+offset) = pftcon%slatop(patch%itype(i)) + ! clm_statevec(cc+2*clm_varsize+offset) = params_inst%kmax(patch%itype(i),1) + ! clm_statevec(cc+3*clm_varsize+offset) = params_inst%kmax(patch%itype(i),2) + ! cc = cc + 1 + ! end do + ! endif + !hcp LAI if(clmupdate_T==1) then error stop "Not implemented: clmupdate_T.eq.1" @@ -540,6 +665,10 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") use clm_time_manager , only : update_DA_nstep use shr_kind_mod , only : r8 => shr_kind_r8 use clm_instMod, only : waterstate_inst + use clm_instMod, only : bgc_vegetation_inst + use pftconMod , only : pftcon + use PatchType , only : patch + use PhotosynthesisMod, only : params_inst implicit none @@ -548,16 +677,28 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2) real(r8), pointer :: h2osoi_ice(:,:) + real(r8), pointer :: leafc(:) ! leaf carbon of patch + real(r8), pointer :: leafn(:) ! leaf nitrogen of patch + +! real(r8) :: tlai(clm_begp:clm_endp) + real(r8) :: incr_lai integer :: i + integer :: cc + integer :: offset character (len = 31) :: fn !TSMP-PDAF: function name for state vector output + character (len = 31) :: fn2 !TSMP-PDAF: function name for state vector outpu character (len = 32) :: fn5 !TSMP-PDAF: function name for state vector outpu character (len = 32) :: fn6 !TSMP-PDAF: function name for state vector outpu + character (len = 32) :: fn7 !TSMP-PDAF: function name for state vector outpu logical :: swc_zero_before_update swc_zero_before_update = .false. + cc = 1 + offset = 0 + #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN ! TSMP-PDAF: For debug runs, output the state vector in files @@ -573,6 +714,11 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") h2osoi_liq => waterstate_inst%h2osoi_liq_col h2osoi_ice => waterstate_inst%h2osoi_ice_col + leafc => bgc_vegetation_inst%cnveg_carbonstate_inst%leafc_patch + leafn => bgc_vegetation_inst%cnveg_nitrogenstate_inst%leafn_patch + + + #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN @@ -619,6 +765,107 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") call update_clm_texture(tstartcycle, mype) endif + ! LAI assimilation: + ! Case 1: + if(clmupdate_lai==1) then + ! all patches + do i=clm_begp,clm_endp + ! DEBUG ONLY + print *, "LST DEBUG leafc,tlai(p) BEFORE:",i,";",leafc(i),";",tlai(i) + + ! At this time tlai(patches) contains the pre-DA patch level LAI + ! clm_statevec(gridcells) contains the post-DA grid level LAI + ! old_frac(patches) contains the fraction of the pre-DA patch lai / grid cell average lai + ! To update tlai(patches) = new gridcell average * old_frac(patches) / patchweight. + ! Then use that tlai to determine leafc(patches) + + !DEBUG + print *, "LST DEBUG old_frac BEFORE", i, ";", old_frac(i) + + if (patch%wtgcell(i) > 0._r8) then ! divide by zero protection + tlai(i) = clm_statevec(patch%gridcell(i)) * old_frac(i)/patch%wtgcell(i) + else + tlai(i) = 0._r8 + endif + ! DEBUG + print *, "LST DEBUG calc tlai BEFORE", i, ";", & + (pftcon%slatop(patch%itype(i))*exp(leafc(i)*pftcon%dsladlai(patch%itype(i)) - 1._r8)) & + /pftcon%dsladlai(patch%itype(i)) + ! update leafc based on Eq 3 from Thornton and Zimmerman, 2007, J Clim, 20, 3902-3923 + ! reformulate the equation to solve for leafc + ! => leafc(p) = log(((tlai(p) * dsladlai(ivt(p)))/slatop(ivt(p))) + 1._r8) / dsladlai(ivt(p)) + if (tlai(i) > 0._r8) then ! invalid log protection + if (pftcon%dsladlai(patch%itype(i)) > 0._r8) then + leafc(i) = log(((tlai(i) * pftcon%dsladlai(patch%itype(i))) / pftcon%slatop(patch%itype(i))) + 1.0_r8) & + / pftcon%dsladlai(patch%itype(i)) + else + leafc(i) = tlai(i) * pftcon%slatop(patch%itype(i)) + endif + else + leafc(i) = 0._r8 + endif + leafc(i) = max(0._r8, leafc(i)) ! Don't allow negative leaf carbon + leafn(i) = leafc(i) / pftcon%leafcn(patch%itype(i)) + + ! DEBUG ONLY + print *, "LST DEBUG leafc,tlai(wt) AFTER: ",i,";",leafc(i),";",tlai(i) + end do + +#ifdef PDAF_DEBUG + IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn7, "(a,i5.5,a,i5.5,a)") "leafc_", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn2, action="write") + WRITE (71,"(es22.15)") leafc(:) + CLOSE(71) + END IF +#endif + endif + + ! Case 2: statevec contains leafc, slatop, dsladlai + ! update of diagnostic tlai left to be done by clm5 itself. + if(clmupdate_lai==2) then + cc = 1 + do i=clm_begp,clm_endp + leafc(i) = clm_statevec(cc) + leafc(i) = max(0._r8, leafc(i)) + leafn(i) = leafc(i) / pftcon%leafcn(patch%itype(i)) + if(clmupdate_lai_params==2) then + pftcon%slatop(patch%itype(i)) = clm_statevec(cc+1*clm_varsize) + pftcon%dsladlai(patch%itype(i)) = clm_statevec(cc+2*clm_varsize) + endif + + cc = cc + 1 + enddo + endif + + if (clmupdate_lai_params==1) then + cc = 1 + do i=clm_begp,clm_endp + pftcon%slatop(patch%itype(i)) = clm_statevec(cc+1*clm_varsize+offset) + ! DEBUG ONLY + print *, "LST DEBUG PARAM,i: ", pftcon%slatop(patch%itype(i)), ",", i + cc = cc + 1 + end do + endif + + ! if (clmupdate_lai_params==3) then + ! cc = 1 + ! do i=clm_begp,clm_endp + ! pftcon%slatop(patch%itype(i)) = clm_statevec(cc+1*clm_varsize+offset) + ! params_inst%kmax(patch%itype(i),1) = clm_statevec(cc+2*clm_varsize+offset) + ! params_inst%kmax(patch%itype(i),2) = clm_statevec(cc+3*clm_varsize+offset) + + ! ! DEBUG ONLY + ! print *, "LST DEBUG PARAM slatop,i: ", pftcon%slatop(patch%itype(i)), ",", i + ! print *, "LST DEBUG PARAM kmax sun,i: ", params_inst%kmax(patch%itype(i),1), ",", i + ! print *, "LST DEBUG PARAM kmax shade,i: ", params_inst%kmax(patch%itype(i),2), ",", i + + ! cc = cc + 1 + ! end do + + ! end if + end subroutine update_clm diff --git a/interface/model/wrapper_tsmp.c b/interface/model/wrapper_tsmp.c index 56687e03..c20d9a4b 100644 --- a/interface/model/wrapper_tsmp.c +++ b/interface/model/wrapper_tsmp.c @@ -198,7 +198,7 @@ void integrate_tsmp() { void update_tsmp(){ #if defined CLMSA - if((model == tag_model_clm) && ((clmupdate_swc != 0) || (clmupdate_T != 0))){ + if((model == tag_model_clm) && ((clmupdate_swc != 0) || (clmupdate_T != 0) || (clmupdate_lai != 0))){ update_clm(&tstartcycle, &mype_world); if(clmprint_swc == 1 || clmupdate_texture == 1 || clmupdate_texture == 2){ print_update_clm(&tcycle, &total_steps);