46subroutine gaussian_to_fv3_interp(lsoil_incr,rla,rlo, &
 
   47                        stcinc,slcinc,soilsnow_tile,lensfc,lsoil,idim,jdim,lsm, myrank)
 
   57    integer, 
intent(in)      :: lsoil_incr, lensfc, lsoil, idim, jdim, myrank, lsm
 
   59    integer, 
intent(in)      :: soilsnow_tile(lensfc)
 
   60    real, 
intent(inout)      :: rla(lensfc), rlo(lensfc)
 
   62    real, 
intent(out)        :: stcinc(lensfc,lsoil)
 
   63    real, 
intent(out)        :: slcinc(lensfc,lsoil)
 
   65    integer                  :: iopt, nret, kgds_gaus(200)
 
   66    integer                  :: igaus, jgaus, ij
 
   68    integer                  :: itile, jtile
 
   70    integer                  :: igausp1, jgausp1
 
   71    logical                  :: upd_slc, upd_stc
 
   74    integer, 
allocatable     :: id1(:,:), id2(:,:), jdc(:,:)
 
   77    real                     :: stc_inc(lsoil)
 
   78    real                     :: slc_inc(lsoil)
 
   79    real, 
allocatable        :: xpts(:), ypts(:), lats(:), lons(:)
 
   80    real, 
allocatable        :: dum2d(:,:), lats_rad(:), lons_rad(:)
 
   81    real, 
allocatable        :: agrid(:,:,:), s2c(:,:,:)
 
   84    logical                  :: gaus_has_soil 
 
   96    kgds_gaus(8)  = nint(-360000./float(
idim_gaus))  
 
   97    kgds_gaus(9)  = nint((360.0 / float(
idim_gaus))*1000.0)
 
  104    if (lsm==lsm_noah) 
then  
  107    elseif (lsm==lsm_noahmp) 
then  
  112    print*,
' Start interpolating first ', lsoil_incr, 
' surface layers only' 
  132    print*,
'fatal error: problem in gdswzd. stop.' 
  133    call mpi_abort(mpi_comm_world, 12, ierr)
 
  136    deallocate (xpts, ypts)
 
  145    lats_rad(j) = dum2d(1,
jdim_gaus-j+1) * 3.1415926 / 180.0
 
  151    lons_rad = dum2d(:,1) * 3.1415926 / 180.0
 
  155    allocate(agrid(idim,jdim,2))
 
  156    agrid(:,:,1) = reshape(rlo, (/idim,jdim/) )
 
  157    agrid(:,:,2) = reshape(rla, (/idim,jdim/) )
 
  158    agrid        = agrid * 3.1415926 / 180.0
 
  160    allocate(id1(idim,jdim))
 
  161    allocate(id2(idim,jdim))
 
  162    allocate(jdc(idim,jdim))
 
  163    allocate(s2c(idim,jdim,4))
 
  172                  lons_rad, lats_rad, id1, id2, jdc, s2c, agrid )
 
  174    deallocate(lons_rad, lats_rad, agrid)
 
  180    ij_loop : 
do ij = 1, lensfc
 
  182        mask_tile    = soilsnow_tile(ij)
 
  186        jtile = (ij-1) / idim + 1
 
  188        if (itile==0) itile = idim
 
  194        if (mask_tile == 1) 
then 
  196           igaus   = id1(itile,jtile)
 
  197           jgaus   = jdc(itile,jtile)
 
  198           igausp1 = id2(itile,jtile)
 
  199           jgausp1 = jdc(itile,jtile)+1
 
  202           gaus_has_soil = .false.
 
  208           if (.not. gaus_has_soil) 
then 
  219                   stc_inc(k)  = stc_inc(k) + (s2c(itile,jtile,1) * 
stc_inc_gaus(k,igaus,jgaus))
 
  221                   slc_inc(k)  = slc_inc(k) + (s2c(itile,jtile,1) * 
slc_inc_gaus(k,igaus,jgaus))
 
  223             wsum  = wsum + s2c(itile,jtile,1)
 
  229                   stc_inc(k) = stc_inc(k) + (s2c(itile,jtile,2) * 
stc_inc_gaus(k,igausp1,jgaus))
 
  231                   slc_inc(k) = slc_inc(k) + (s2c(itile,jtile,2) * 
slc_inc_gaus(k,igausp1,jgaus))
 
  233             wsum  = wsum + s2c(itile,jtile,2)
 
  239                   stc_inc(k) = stc_inc(k) + (s2c(itile,jtile,3) * 
stc_inc_gaus(k,igausp1,jgausp1))
 
  241                   slc_inc(k) = slc_inc(k) + (s2c(itile,jtile,3) * 
slc_inc_gaus(k,igausp1,jgausp1))
 
  243             wsum  = wsum + s2c(itile,jtile,3)
 
  249                   stc_inc(k) = stc_inc(k) + (s2c(itile,jtile,4) * 
stc_inc_gaus(k,igaus,jgausp1))
 
  251                   slc_inc(k) = slc_inc(k) + (s2c(itile,jtile,4) * 
slc_inc_gaus(k,igaus,jgausp1))
 
  253             wsum  = wsum + s2c(itile,jtile,4)
 
  258             stc_inc(k) = stc_inc(k) / wsum
 
  259             stcinc(ij,k) = stc_inc(k)
 
  260             slc_inc(k) = slc_inc(k) / wsum
 
  261             slcinc(ij,k) = slc_inc(k)
 
  268    write(*,
'(a,i2)') 
' Finish soil increments interpolation for rank : ', myrank
 
  270    deallocate(id1, id2, jdc, s2c)
 
 
  295subroutine add_increment_soil(lsoil_incr,stcinc,slcinc,stc_state,smc_state,slc_state,stc_updated,&
 
  296              slc_updated,soilsnow_tile,soilsnow_fg_tile,lensfc,lsoil,lsm,myrank)
 
  302    integer, 
intent(in)      :: lsoil_incr, lensfc, lsoil, myrank, lsm
 
  304    integer, 
intent(in)      :: soilsnow_tile(lensfc), soilsnow_fg_tile(lensfc)
 
  305    real, 
intent(inout)      :: stc_state(lensfc, lsoil)
 
  306    real, 
intent(inout)      :: slc_state(lensfc, lsoil)
 
  307    real, 
intent(inout)      :: smc_state(lensfc, lsoil)
 
  308    integer, 
intent(out)     :: stc_updated(lensfc), slc_updated(lensfc)
 
  312    integer                  :: mask_tile, mask_fg_tile
 
  313    logical                  :: upd_slc, upd_stc
 
  315    real                     :: stcinc(lensfc,lsoil)
 
  316    real                     :: slcinc(lensfc,lsoil)
 
  318    integer                  :: k, nother, nsnowupd
 
  319    integer                  :: nstcupd, nslcupd,  nfrozen, nfrozen_upd
 
  320    logical                  :: soil_freeze, soil_ice
 
  325    if (lsm==lsm_noah) 
then 
  328    elseif (lsm==lsm_noahmp) 
then 
  334    print*,
'adjust soil using increments on cubed-sphere tiles' 
  335    print*,
'updating soil temps', upd_stc
 
  336    print*,
'updating soil moisture', upd_slc
 
  337    print*,
'adjusting first ', lsoil_incr, 
' surface layers only' 
  347    ij_loop : 
do ij = 1, lensfc
 
  349        mask_tile    = soilsnow_tile(ij)
 
  350        mask_fg_tile = soilsnow_fg_tile(ij)
 
  356        if (mask_tile <= 0) 
then  
  365        if (mask_fg_tile == 2 .or. mask_tile == 2) 
then 
  366         nsnowupd = nsnowupd + 1
 
  374        if (mask_tile == 1) 
then 
  384             if ( stc_state(ij,k) < tfreez)  soil_freeze=.true.
 
  385             if ( smc_state(ij,k) - slc_state(ij,k) > 0.001 )  soil_ice=.true.
 
  388                stc_state(ij,k) = stc_state(ij,k) + stcinc(ij,k)
 
  391                    nstcupd = nstcupd + 1
 
  395             if ( (stc_state(ij,k) < tfreez) .and. (.not. soil_freeze) .and. (k==1) )&
 
  396                   nfrozen_upd = nfrozen_upd + 1
 
  399             if ( (.not. soil_freeze ) .and. (.not. soil_ice ) ) 
then 
  402                    nslcupd = nslcupd + 1
 
  407                   slc_state(ij,k) = max(slc_state(ij,k) + slcinc(ij,k), 0.0)
 
  408                   smc_state(ij,k) = max(smc_state(ij,k) + slcinc(ij,k), 0.0)
 
  411                if (k==1) nfrozen = nfrozen+1
 
  420   write(*,
'(a,i2)') 
' statistics of grids number processed for rank : ', myrank
 
  421   write(*,
'(a,i8)') 
' soil grid total', lensfc
 
  422   write(*,
'(a,i8)') 
' soil grid cells slc updated = ',nslcupd
 
  423   write(*,
'(a,i8)') 
' soil grid cells stc updated = ',nstcupd
 
  424   write(*,
'(a,i8)') 
' soil grid cells not updated, frozen = ',nfrozen
 
  425   write(*,
'(a,i8)') 
' soil grid cells update, became frozen = ',nfrozen_upd
 
  426   write(*,
'(a,i8)') 
' (not updated yet) snow grid cells = ', nsnowupd
 
  427   write(*,
'(a,i8)') 
' grid cells, without soil or snow = ', nother
 
 
  533subroutine apply_land_da_adjustments_soil(lsoil_incr, lsm, isot, ivegsrc,lensfc, &
 
  534                 lsoil, rsoiltype, mask, stc_bck, stc_adj, smc_adj, slc_adj, &
 
  535                 stc_updated, slc_updated, zsoil)
 
  538    use set_soilveg_snippet_mod, 
only: set_soilveg_noah,set_soilveg_noahmp
 
  539    use sflx_snippet,    
only: frh2o
 
  543    integer, 
intent(in)           :: lsoil_incr, lsm, lensfc, lsoil, isot, ivegsrc
 
  544    real, 
intent(in)              :: rsoiltype(lensfc) 
 
  545    integer, 
intent(in)           :: mask(lensfc)
 
  546    real, 
intent(in)              :: stc_bck(lensfc, lsoil)
 
  547    integer, 
intent(in)           :: stc_updated(lensfc), slc_updated(lensfc)
 
  548    real, 
intent(inout)           :: smc_adj(lensfc,lsoil), slc_adj(lensfc,lsoil) 
 
  549    real, 
intent(inout)           :: stc_adj(lensfc, lsoil)
 
  550    real(kind=4), intent(in)      :: zsoil(lsoil)
 
  553    logical                       :: frzn_bck, frzn_anl
 
  554    logical                       :: soil_freeze, soil_ice
 
  556    integer                       :: i, l, n_freeze, n_thaw, ierr
 
  557    integer                       :: myrank, soiltype, iret, n_stc, n_slc
 
  558    logical                       :: upd_slc, upd_stc
 
  562    real, 
parameter               :: tfreez=273.16
 
  563    real, 
dimension(30)           :: maxsmc, bb, satpsi
 
  564    real, 
dimension(4)            :: dz 
 
  566    real, 
parameter          :: hfus=0.3336e06 
 
  567    real, 
parameter          :: grav=9.80616
 
  570    call mpi_comm_rank(mpi_comm_world, myrank, ierr)
 
  572    if (lsm==lsm_noah) 
then  
  575    elseif (lsm==lsm_noahmp) 
then  
  583        call set_soilveg_noah(isot, ivegsrc, maxsmc, bb, satpsi, iret)
 
  585            print *, 
'FATAL ERROR: problem in set_soilveg_noah' 
  586            call mpi_abort(mpi_comm_world, 10, ierr)
 
  589        print *, 
'Adjusting noah model smc after stc DA update' 
  595          if (mask(i) > 0) 
then  
  597               frzn_bck = (stc_bck(i,l) .LT. tfreez )
 
  598               frzn_anl = (stc_adj(i,l) .LT. tfreez )
 
  600               if (frzn_bck .eqv. frzn_anl) 
then 
  602               elseif (frzn_bck .and. .not. frzn_anl) 
then 
  605                    n_freeze = n_freeze + 1
 
  609               soiltype = nint(rsoiltype(i))
 
  611               call frh2o(stc_adj(i,l), smc_adj(i,l),slc_adj(i,l), maxsmc(soiltype), &
 
  612                          bb(soiltype), satpsi(soiltype),slc_new)
 
  614               slc_adj(i,l) = max( min( slc_new, smc_adj(i,l)), 0.0 )
 
  618        print *, 
'adjusted: ', n_thaw,
' thawed,', n_freeze, 
' frozen' 
  624          call set_soilveg_noahmp(isot, ivegsrc, maxsmc, bb, satpsi, iret)
 
  626               print *, 
'FATAL ERROR: problem in set_soilveg_noahmp' 
  627               call mpi_abort(mpi_comm_world, 10, ierr)
 
  634            if (stc_updated(i) == 1 ) 
then  
  636                soiltype = nint(rsoiltype(i))
 
  641                   if (stc_adj(i,l) .LT. tfreez )
then 
  643                      smp = hfus*(tfreez-stc_adj(i,l))/(grav*stc_adj(i,l)) 
 
  644                      slc_new=maxsmc(soiltype)*(smp/satpsi(soiltype))**(-1./bb(soiltype))
 
  645                      slc_adj(i,l) = max( min( slc_new, smc_adj(i,l)), 0.0 )
 
  648                   if (stc_adj(i,l) .GT. tfreez )
then  
  649                      slc_adj(i,l)=smc_adj(i,l)
 
  661              dz(l) = -zsoil(l) + zsoil(l-1) 
 
  663          print *, 
'Applying soil moisture mins '  
  666          if (slc_updated(i) == 1 ) 
then  
  672                slc_adj(i,l) = max( 0.001/dz(l), slc_adj(i,l) )
 
  673                smc_adj(i,l) = max( 0.001/dz(l), smc_adj(i,l) )
 
  680        print *, 
'FATAL ERROR: unrecognised LSM,', lsm
 
  681        call mpi_abort(mpi_comm_world, 10, ierr)
 
  684    write(*,
'(a,i2)') 
'statistics of grids number processed for rank : ', myrank 
 
  685    write(*,
'(a,i8)') 
' soil grid total', lensfc
 
  686    write(*,
'(a,i8)') 
' soil grid cells with slc update', n_slc
 
  687    write(*,
'(a,i8)') 
' soil grid cells with stc update', n_stc
 
 
  706subroutine apply_land_da_adjustments_snd(lsm, lensfc, mask, swe_bck, snd_bck, snd_anl, swe_adj)
 
  709    use bulk_snow_module, 
only: calc_density
 
  713    integer, 
intent(in) :: lsm, lensfc
 
  714    integer, 
intent(in) :: mask(lensfc)
 
  715    real, 
intent(in)    :: swe_bck(lensfc), snd_bck(lensfc)
 
  716    real, 
intent(in)    :: snd_anl(lensfc)
 
  717    real, 
intent(inout)    :: swe_adj(lensfc)
 
  719    integer  :: ierr, myrank, i
 
  721    real     :: density(lensfc)
 
  723    call mpi_comm_rank(mpi_comm_world, myrank, ierr)
 
  725    if (lsm .NE. lsm_noah) 
then 
  726        print *, 
'FATAL ERROR: apply_land_da_adjustments not coded for models other than noah', lsm
 
  727        call mpi_abort(mpi_comm_world, 10, ierr)
 
  731   call calc_density(lensfc, mask, swe_bck, snd_bck, myrank, density)
 
  734        if ( mask(i)>0 ) 
then 
  735                swe_adj(i) = snd_anl(i)*density(i)
 
 
subroutine, public remap_coef(is, ie, js, je, im, jm, lon, lat, id1, id2, jdc, s2c, agrid)
Generate the weights and index of the grids used in the bilinear interpolation.