global_cycle  1.11.0
land_increments.f90
Go to the documentation of this file.
1 
4 
5 module land_increments
6 
7  private
8 
9  public add_increment_soil
10  public add_increment_snow
11  public calculate_landinc_mask
12  public apply_land_da_adjustments_soil
13  public apply_land_da_adjustments_snd
14  public lsm_noah, lsm_noahmp
15 
16  integer, parameter :: lsm_noah=1
17  integer, parameter :: lsm_noahmp=2
19 
20  ! control state for soil analysis:
21  integer, parameter :: lsoil_incr=3
22 
23  real, parameter :: tfreez=273.16
24 contains
25 
54 
55 subroutine add_increment_soil(rla,rlo,stc_state,smc_state,slc_state,stc_updated, slc_updated, &
56  soilsnow_tile,soilsnow_fg_tile,lensfc,lsoil,idim,jdim,lsm, myrank)
57 
58  use utils
59  use gdswzd_mod
60  use read_write_data, only : idim_gaus, jdim_gaus, &
61  stc_inc_gaus, soilsnow_gaus, slc_inc_gaus
62  use mpi
63 
64  implicit none
65 
66  integer, intent(in) :: lensfc, lsoil, idim, jdim, myrank, lsm
67 
68  integer, intent(in) :: soilsnow_tile(lensfc), soilsnow_fg_tile(lensfc)
69  real, intent(inout) :: rla(lensfc), rlo(lensfc)
70  real, intent(inout) :: stc_state(lensfc, lsoil)
71  real, intent(inout) :: slc_state(lensfc, lsoil)
72  real, intent(inout) :: smc_state(lensfc, lsoil)
73  integer, intent(out) :: stc_updated(lensfc), slc_updated(lensfc)
74 
75  integer :: iopt, nret, kgds_gaus(200)
76  integer :: igaus, jgaus, ij
77  integer :: mask_tile, mask_fg_tile
78  integer :: itile, jtile
79  integer :: j, ierr
80  integer :: igausp1, jgausp1
81  logical :: upd_slc, upd_stc
82  real :: fill
83 
84  integer, allocatable :: id1(:,:), id2(:,:), jdc(:,:)
85 
86  real :: wsum
87  real :: stc_inc(lsoil)
88  real :: slc_inc(lsoil)
89  real, allocatable :: xpts(:), ypts(:), lats(:), lons(:)
90  real, allocatable :: dum2d(:,:), lats_rad(:), lons_rad(:)
91  real, allocatable :: agrid(:,:,:), s2c(:,:,:)
92 
93  integer :: k, nother, nsnowupd, nnosoilnear
94  integer :: nstcupd, nslcupd, nfrozen, nfrozen_upd
95  logical :: gaus_has_soil, soil_freeze, soil_ice
96 
97  stc_updated=0
98  slc_updated=0
99  ! this produces the same lat/lon as can be read in from the file
100 
101  kgds_gaus = 0
102  kgds_gaus(1) = 4 ! oct 6 - type of grid (gaussian)
103  kgds_gaus(2) = idim_gaus ! oct 7-8 - # pts on latitude circle
104  kgds_gaus(3) = jdim_gaus
105  kgds_gaus(4) = 90000 ! oct 11-13 - lat of origin
106  kgds_gaus(5) = 0 ! oct 14-16 - lon of origin
107  kgds_gaus(6) = 128 ! oct 17 - resolution flag
108  kgds_gaus(7) = -90000 ! oct 18-20 - lat of extreme point
109  kgds_gaus(8) = nint(-360000./float(idim_gaus)) ! oct 21-23 - lon of extreme point
110  kgds_gaus(9) = nint((360.0 / float(idim_gaus))*1000.0)
111  ! oct 24-25 - longitude direction incr.
112  kgds_gaus(10) = jdim_gaus/2 ! oct 26-27 - number of circles pole to equator
113  kgds_gaus(12) = 255 ! oct 29 - reserved
114  kgds_gaus(20) = 255 ! oct 5 - not used, set to 255
115 
116 
117  if (lsm==lsm_noah) then
118  upd_stc=.true.
119  upd_slc=.false. ! not coded
120  elseif (lsm==lsm_noahmp) then
121  upd_stc=.true.
122  upd_slc=.true.
123  endif
124 
125  print*
126  print*,'adjust soil using gsi increments on gaussian grid'
127  print*,'updating soil temps', upd_stc
128  print*,'updating soil moisture', upd_slc
129  print*,'adjusting first ', lsoil_incr, ' surface layers only'
130 
131  !----------------------------------------------------------------------
132  ! call gdswzd to compute the lat/lon of each gsi gaussian grid point.
133  !----------------------------------------------------------------------
134 
135  iopt = 0
136  fill = -9999.
137  allocate(xpts(idim_gaus*jdim_gaus))
138  allocate(ypts(idim_gaus*jdim_gaus))
139  allocate(lats(idim_gaus*jdim_gaus))
140  allocate(lons(idim_gaus*jdim_gaus))
141  xpts = fill
142  ypts = fill
143  lats = fill
144  lons = fill
145 
146  call gdswzd(kgds_gaus,iopt,(idim_gaus*jdim_gaus),fill,xpts,ypts,lons,lats,nret)
147 
148  if (nret /= (idim_gaus*jdim_gaus)) then
149  print*,'fatal error: problem in gdswzd. stop.'
150  call mpi_abort(mpi_comm_world, 12, ierr)
151  endif
152 
153  deallocate (xpts, ypts)
154 
155  allocate(dum2d(idim_gaus,jdim_gaus))
156  dum2d = reshape(lats, (/idim_gaus,jdim_gaus/) )
157  deallocate(lats)
158 
159  allocate(lats_rad(jdim_gaus))
160 
161  do j = 1, jdim_gaus
162  lats_rad(j) = dum2d(1,jdim_gaus-j+1) * 3.1415926 / 180.0
163  enddo
164 
165  dum2d = reshape(lons, (/idim_gaus,jdim_gaus/) )
166  deallocate(lons)
167  allocate(lons_rad(idim_gaus))
168  lons_rad = dum2d(:,1) * 3.1415926 / 180.0
169 
170  deallocate(dum2d)
171 
172  allocate(agrid(idim,jdim,2))
173  agrid(:,:,1) = reshape(rlo, (/idim,jdim/) )
174  agrid(:,:,2) = reshape(rla, (/idim,jdim/) )
175  agrid = agrid * 3.1415926 / 180.0
176 
177  allocate(id1(idim,jdim))
178  allocate(id2(idim,jdim))
179  allocate(jdc(idim,jdim))
180  allocate(s2c(idim,jdim,4))
181 
182  !----------------------------------------------------------------------
183  ! compute bilinear weights for each model point from the nearest
184  ! four gsi/gaussian points. does not account for mask. that
185  ! happens later.
186  !----------------------------------------------------------------------
187 
188  call remap_coef( 1, idim, 1, jdim, idim_gaus, jdim_gaus, &
189  lons_rad, lats_rad, id1, id2, jdc, s2c, agrid )
190 
191  deallocate(lons_rad, lats_rad, agrid)
192  !
193  ! initialize variables for counts statitics to be zeros
194  !
195 
196  !
197  nother = 0 ! grid cells not land
198  nsnowupd = 0 ! grid cells with snow (temperature not yet updated)
199  nnosoilnear = 0 ! grid cells where model has soil, but 4 closest gaus grids don't
200  ! (no update made here)
201  nslcupd = 0 ! grid cells that are updated
202  nstcupd = 0 ! grid cells that are updated
203  nfrozen = 0 ! not update as frozen soil
204  nfrozen_upd = 0 ! not update as frozen soil
205 
206 
207  ij_loop : do ij = 1, lensfc
208 
209  mask_tile = soilsnow_tile(ij)
210  mask_fg_tile = soilsnow_fg_tile(ij)
211 
212  !----------------------------------------------------------------------
213  ! mask: 1 - soil, 2 - snow, 0 - land-ice, -1 - not land
214  !----------------------------------------------------------------------
215 
216  if (mask_tile <= 0) then ! skip if neither soil nor snow
217  nother = nother + 1
218  cycle ij_loop
219  endif
220 
221 
222  ! get i,j index on array of (idim,jdim) from known ij
223 
224  jtile = (ij-1) / idim + 1
225  itile = mod(ij,idim)
226  if (itile==0) itile = idim
227 
228  !----------------------------------------------------------------------
229  ! if snow is present before or after snow update, skip soil analysis
230  !----------------------------------------------------------------------
231 
232  if (mask_fg_tile == 2 .or. mask_tile == 2) then
233  nsnowupd = nsnowupd + 1
234  cycle ij_loop
235  endif
236 
237  !----------------------------------------------------------------------
238  ! do update to soil temperature grid cells, using bilinear interp
239  !----------------------------------------------------------------------
240 
241  if (mask_tile == 1) then
242  ! these are the four nearest grid cells on the gaus grid
243  igaus = id1(itile,jtile)
244  jgaus = jdc(itile,jtile)
245  igausp1 = id2(itile,jtile)
246  jgausp1 = jdc(itile,jtile)+1
247 
248  ! make sure gaus grid has soil nearby
249  gaus_has_soil = .false.
250  if (soilsnow_gaus(igaus,jgaus) == 1 .or. &
251  soilsnow_gaus(igausp1,jgaus) == 1 .or. &
252  soilsnow_gaus(igausp1,jgausp1) == 1 .or. &
253  soilsnow_gaus(igaus,jgausp1) == 1) gaus_has_soil = .true.
254 
255  if (.not. gaus_has_soil) then
256  nnosoilnear = nnosoilnear + 1
257  cycle ij_loop
258  endif
259 
260  stc_inc = 0.0
261  slc_inc = 0.0
262  wsum = 0.0
263 
264  if (soilsnow_gaus(igaus,jgaus) == 1) then
265  do k = 1, lsoil_incr
266  if (upd_stc) &
267  stc_inc(k) = stc_inc(k) + (s2c(itile,jtile,1) * stc_inc_gaus(k,igaus,jgaus))
268  if (upd_slc) &
269  slc_inc(k) = slc_inc(k) + (s2c(itile,jtile,1) * slc_inc_gaus(k,igaus,jgaus))
270  enddo
271  wsum = wsum + s2c(itile,jtile,1)
272  endif
273 
274  if (soilsnow_gaus(igausp1,jgaus) == 1) then
275  do k = 1, lsoil_incr
276  if (upd_stc) &
277  stc_inc(k) = stc_inc(k) + (s2c(itile,jtile,2) * stc_inc_gaus(k,igausp1,jgaus))
278  if (upd_slc) &
279  slc_inc(k) = slc_inc(k) + (s2c(itile,jtile,2) * slc_inc_gaus(k,igausp1,jgaus))
280  enddo
281  wsum = wsum + s2c(itile,jtile,2)
282  endif
283 
284  if (soilsnow_gaus(igausp1,jgausp1) == 1) then
285  do k = 1, lsoil_incr
286  if (upd_stc) &
287  stc_inc(k) = stc_inc(k) + (s2c(itile,jtile,3) * stc_inc_gaus(k,igausp1,jgausp1))
288  if (upd_slc) &
289  slc_inc(k) = slc_inc(k) + (s2c(itile,jtile,3) * slc_inc_gaus(k,igausp1,jgausp1))
290  enddo
291  wsum = wsum + s2c(itile,jtile,3)
292  endif
293 
294  if (soilsnow_gaus(igaus,jgausp1) == 1) then
295  do k = 1, lsoil_incr
296  if (upd_stc) &
297  stc_inc(k) = stc_inc(k) + (s2c(itile,jtile,4) * stc_inc_gaus(k,igaus,jgausp1))
298  if (upd_slc) &
299  slc_inc(k) = slc_inc(k) + (s2c(itile,jtile,4) * slc_inc_gaus(k,igaus,jgausp1))
300  enddo
301  wsum = wsum + s2c(itile,jtile,4)
302  endif
303 
304  ! normalize increments
305  do k = 1, lsoil_incr
306  stc_inc(k) = stc_inc(k) / wsum
307  slc_inc(k) = slc_inc(k) / wsum
308  enddo
309  !----------------------------------------------------------------------
310  ! add the interpolated increment to the background
311  !----------------------------------------------------------------------
312 
313  soil_freeze=.false.
314  soil_ice=.false.
315  do k = 1, lsoil_incr
316 
317  if ( stc_state(ij,k) < tfreez) soil_freeze=.true.
318  if ( smc_state(ij,k) - slc_state(ij,k) > 0.001 ) soil_ice=.true.
319 
320  if (upd_stc) then
321  stc_state(ij,k) = stc_state(ij,k) + stc_inc(k)
322  if (k==1) then
323  stc_updated(ij) = 1
324  nstcupd = nstcupd + 1
325  endif
326  endif
327 
328  if ( (stc_state(ij,k) < tfreez) .and. (.not. soil_freeze) .and. (k==1) ) &
329  nfrozen_upd = nfrozen_upd + 1
330 
331  ! do not do updates if this layer or any above is frozen
332  if ( (.not. soil_freeze ) .and. (.not. soil_ice ) ) then
333  if (upd_slc) then
334  if (k==1) then
335  nslcupd = nslcupd + 1
336  slc_updated(ij) = 1
337  endif
338  ! apply zero limit here (higher, model-specific limits are later)
339  slc_state(ij,k) = max(slc_state(ij,k) + slc_inc(k), 0.0)
340  smc_state(ij,k) = max(smc_state(ij,k) + slc_inc(k), 0.0)
341  endif
342  else
343  if (k==1) nfrozen = nfrozen+1
344  endif
345 
346  enddo
347 
348  endif ! if soil/snow point
349 
350  enddo ij_loop
351 
352  write(*,'(a,i2)') 'statistics of grids number processed for rank : ', myrank
353  write(*,'(a,i8)') ' soil grid total', lensfc
354  write(*,'(a,i8)') ' soil grid cells slc updated = ',nslcupd
355  write(*,'(a,i8)') ' soil grid cells stc updated = ',nstcupd
356  write(*,'(a,i8)') ' soil grid cells not updated, frozen = ',nfrozen
357  write(*,'(a,i8)') ' soil grid cells update, became frozen = ',nfrozen_upd
358  write(*,'(a,i8)') ' (not updated) soil grid cells, no soil nearby on gsi grid = ',nnosoilnear
359  write(*,'(a,i8)') ' (not updated yet) snow grid cells = ', nsnowupd
360  write(*,'(a,i8)') ' grid cells, without soil or snow = ', nother
361 
362  deallocate(id1, id2, jdc, s2c)
363 
364 end subroutine add_increment_soil
365 
377 
378 subroutine add_increment_snow(snd_inc,mask,lensfc,snd)
379 
380  implicit none
381 
382  integer, intent(in) :: lensfc
383  real, intent(in) :: snd_inc(lensfc)
384  integer, intent(in) :: mask(lensfc)
385  real, intent(inout) :: snd(lensfc)
386 
387  integer :: i
388 
389 
390  do i =1, lensfc
391  if (mask(i) > 0) then ! if land
392  snd(i) = max( snd(i) + snd_inc(i), 0.)
393  endif
394  enddo
395 
396 end subroutine add_increment_snow
397 
408 subroutine calculate_landinc_mask(smc,swe,vtype,lensfc,veg_type_landice,mask)
409 
410  implicit none
411 
412  integer, intent(in) :: lensfc, veg_type_landice
413  real, intent(in) :: smc(lensfc), swe(lensfc)
414  real, intent(in) :: vtype(lensfc)
415  integer, intent(out) :: mask(lensfc)
416 
417  integer :: i
418 
419  mask = -1 ! not land
420 
421  ! land (but not land-ice)
422  do i=1,lensfc
423  if (smc(i) .LT. 0.99) then
424  if (swe(i) .GT. 0.001) then ! snow covered land
425  mask(i) = 2
426  else ! non-snow covered land
427  mask(i) = 1
428  endif
429  end if ! else should work here too
430  if ( nint(vtype(i)) == veg_type_landice ) then ! land-ice
431  mask(i) = 0
432  endif
433  end do
434 
435 end subroutine calculate_landinc_mask
436 
464 
465 subroutine apply_land_da_adjustments_soil(lsm, isot, ivegsrc,lensfc, &
466  lsoil, rsoiltype, mask, stc_bck, stc_adj, smc_adj, slc_adj, &
467  stc_updated, slc_updated, zsoil)
468 
469  use mpi
470  use set_soilveg_snippet_mod, only: set_soilveg
471  use sflx_snippet, only: frh2o
472 
473  implicit none
474 
475  integer, intent(in) :: lsm, lensfc, lsoil, isot, ivegsrc
476  real, intent(in) :: rsoiltype(lensfc) ! soil types, as real
477  integer, intent(in) :: mask(lensfc)
478  real, intent(in) :: stc_bck(lensfc, lsoil)
479  integer, intent(in) :: stc_updated(lensfc), slc_updated(lensfc)
480  real, intent(inout) :: smc_adj(lensfc,lsoil), slc_adj(lensfc,lsoil)
481  real, intent(inout) :: stc_adj(lensfc, lsoil)
482  real(kind=4), intent(in) :: zsoil(lsoil)
483 
484 
485  logical :: frzn_bck, frzn_anl
486  logical :: soil_freeze, soil_ice
487 
488  integer :: i, l, n_freeze, n_thaw, ierr, n_revert
489  integer :: myrank, soiltype, iret, n_stc, n_slc
490  logical :: upd_slc, upd_stc
491 
492  real :: slc_new
493 
494  real, parameter :: tfreez=273.16
495  real, dimension(30) :: maxsmc, bb, satpsi
496  real, dimension(4) :: dz ! layer thickness
497 
498  call mpi_comm_rank(mpi_comm_world, myrank, ierr)
499 
500  if (lsm==lsm_noah) then
501  upd_stc=.true.
502  upd_slc=.false.
503  elseif (lsm==lsm_noahmp) then
504  upd_stc=.true.
505  upd_slc=.true.
506  endif
507 
508  select case (lsm )
509  case(lsm_noah)
510  ! initialise soil properties
511  call set_soilveg(isot, ivegsrc, maxsmc, bb, satpsi, iret)
512  if (iret < 0) then
513  print *, 'FATAL ERROR: problem in set_soilveg'
514  call mpi_abort(mpi_comm_world, 10, ierr)
515  endif
516 
517  print *, 'Adjusting noah model smc after stc DA update'
518 
519  n_freeze = 0
520  n_thaw = 0
521 
522  do i=1,lensfc
523  if (mask(i) > 0) then ! if soil location
524  do l = 1, lsoil
525  frzn_bck = (stc_bck(i,l) .LT. tfreez )
526  frzn_anl = (stc_adj(i,l) .LT. tfreez )
527 
528  if (frzn_bck .eqv. frzn_anl) then
529  cycle
530  elseif (frzn_bck .and. .not. frzn_anl) then
531  n_thaw = n_thaw + 1
532  else
533  n_freeze = n_freeze + 1
534  endif
535 
536  ! make adjustment (same routine for both)
537  soiltype = nint(rsoiltype(i))
538  ! bb and maxsmc are in the namelist_soilveg, need soiltype index
539  call frh2o(stc_adj(i,l), smc_adj(i,l),slc_adj(i,l), maxsmc(soiltype), &
540  bb(soiltype), satpsi(soiltype),slc_new)
541 
542  slc_adj(i,l) = max( min( slc_new, smc_adj(i,l)), 0.0 )
543  enddo
544  endif
545  enddo
546  print *, 'adjusted: ', n_thaw,' thawed,', n_freeze, ' frozen'
547 
548  case (lsm_noahmp)
549 
550  if (upd_stc) then
551 
552  print *, 'Reverting frozen noah-mp model stc back to background'
553  n_revert=0
554  n_stc = 0
555  n_slc = 0
556 
557  do i=1,lensfc
558  if (stc_updated(i) == 1 ) then
559  n_stc = n_stc+1
560  ! remove soil temperature increments if frozen
561  soil_freeze=.false.
562  soil_ice=.false.
563  do l = 1, lsoil_incr
564  if ( min(stc_bck(i,l),stc_adj(i,l)) < tfreez) soil_freeze=.true.
565  if ( smc_adj(i,l) - slc_adj(i,l) > 0.001 ) soil_ice=.true.
566  if ( soil_freeze .or. soil_ice ) then
567  ! for now, revert update. Later, adjust SMC/SLC for update.
568  if (l==1) n_revert = n_revert+1
569  stc_adj(i,l)=stc_bck(i,l)
570  endif
571  enddo
572  endif
573  enddo
574 
575  endif
576  if (upd_slc) then
577 
578  dz(1) = -zsoil(1)
579  do l = 2,lsoil
580  dz(l) = -zsoil(l) + zsoil(l-1)
581  enddo
582  print *, 'Applying soil moisture mins '
583 
584  do i=1,lensfc
585  if (slc_updated(i) == 1 ) then
586  n_slc = n_slc+1
587  ! apply SM bounds (later: add upper SMC limit)
588  do l = 1, lsoil_incr
589  ! noah-mp minimum is 1 mm per layer (in SMC)
590  ! no need to maintain frozen amount, would be v. small.
591  slc_adj(i,l) = max( 0.001/dz(l), slc_adj(i,l) )
592  smc_adj(i,l) = max( 0.001/dz(l), smc_adj(i,l) )
593  enddo
594  endif
595  enddo
596  endif
597 
598  case default
599  print *, 'FATAL ERROR: unrecognised LSM,', lsm
600  call mpi_abort(mpi_comm_world, 10, ierr)
601  end select
602 
603  write(*,'(a,i2)') 'statistics of grids number processed for rank : ', myrank
604  write(*,'(a,i8)') ' soil grid total', lensfc
605  write(*,'(a,i8)') ' soil grid cells with slc update', n_slc
606  write(*,'(a,i8)') ' soil grid cells with stc update', n_stc
607  write(*,'(a,i8)') ' soil grid cells reverted', n_revert
608 
609 end subroutine apply_land_da_adjustments_soil
610 
616 
625 
626 subroutine apply_land_da_adjustments_snd(lsm, lensfc, mask, swe_bck, snd_bck, snd_anl, swe_adj)
627 
628  use mpi
629  use bulk_snow_module, only: calc_density
630 
631  implicit none
632 
633  integer, intent(in) :: lsm, lensfc
634  integer, intent(in) :: mask(lensfc)
635  real, intent(in) :: swe_bck(lensfc), snd_bck(lensfc)
636  real, intent(in) :: snd_anl(lensfc)
637  real, intent(inout) :: swe_adj(lensfc)
638 
639  integer :: ierr, myrank, i
640 
641  real :: density(lensfc)
642 
643  call mpi_comm_rank(mpi_comm_world, myrank, ierr)
644 
645  if (lsm .NE. lsm_noah) then
646  print *, 'FATAL ERROR: apply_land_da_adjustments not coded for models other than noah', lsm
647  call mpi_abort(mpi_comm_world, 10, ierr)
648  endif
649 
650  ! calculate snow density from forecasts
651  call calc_density(lensfc, mask, swe_bck, snd_bck, myrank, density)
652 
653  do i =1, lensfc
654  if ( mask(i)>0 ) then
655  swe_adj(i) = snd_anl(i)*density(i)
656  endif
657  enddo
658 
659 
660 end subroutine apply_land_da_adjustments_snd
661 
662 end module land_increments
utils
Module containing utility routines.
Definition: utils.F90:7