Skip to content

Commit

Permalink
Changes for cloud control variables with MRI-4DVAR,
Browse files Browse the repository at this point in the history 
 tested for cloud_cv_options=3 and use_cv_w=true with radar data

	modified:   var/da/da_main/da_solve.inc
	modified:   var/da/da_recursive_filter/da_transform_through_rf_inv.inc
	modified:   var/da/da_setup_structures/da_write_vp.inc
	modified:   var/da/da_vtox_transforms/da_transform_vptox_inv.inc
	modified:   var/da/da_vtox_transforms/da_transform_vtovv_inv.inc
	modified:   var/da/da_vtox_transforms/da_transform_vtox_inv.inc
	modified:   var/mri4dvar/da_bilin.f90
	modified:   var/mri4dvar/da_vp_bilin.f90
	modified:   var/mri4dvar/da_vp_split.f90
	modified:   var/mri4dvar/run_mri3d4dvar.csh_pbs
  • Loading branch information
@liujake
liujake committed Jul 24, 2017
1 parent 5e94060 commit 606ac0e
Show file tree
Hide file tree
Showing 10 changed files with 441 additions and 334 deletions.
39 changes: 39 additions & 0 deletions var/da/da_main/da_solve.inc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@

real, allocatable :: grid_box_area(:,:), mapfac(:,:)
real, allocatable :: v1(:,:,:),v2(:,:,:),v3(:,:,:),v4(:,:,:),v5(:,:,:)
real, allocatable :: v6(:,:,:),v7(:,:,:),v8(:,:,:),v9(:,:,:),v10(:,:,:),v11(:,:,:)
character (len=10) :: variable_name

integer :: iwin, num_subtwindow
Expand Down Expand Up @@ -608,22 +609,52 @@
allocate( v3(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1) )
allocate( v4(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1) )
allocate( v5(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1) )
if ( cloud_cv_options >= 2 ) then
allocate( v6(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1) )
allocate( v7(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1) )
allocate( v8(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1) )
allocate( v9(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1) )
allocate( v10(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1) )
end if
if ( use_cv_w ) allocate( v11(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1) )

read(vp_unit) v1, v2, v3, v4, v5
if ( cloud_cv_options >= 2 ) read(vp_unit) v6, v7, v8, v9, v10
if ( use_cv_w ) read(vp_unit) v11

if ( use_interpolate_cvt ) then
grid%vv%v1(ips:ipe,jps:jpe,kps:kpe) = v1(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
grid%vv%v2(ips:ipe,jps:jpe,kps:kpe) = v2(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
grid%vv%v3(ips:ipe,jps:jpe,kps:kpe) = v3(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
grid%vv%v4(ips:ipe,jps:jpe,kps:kpe) = v4(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
grid%vv%v5(ips:ipe,jps:jpe,kps:kpe) = v5(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
if ( cloud_cv_options >= 2 ) then
grid%vv%v6(ips:ipe,jps:jpe,kps:kpe) = v6(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
grid%vv%v7(ips:ipe,jps:jpe,kps:kpe) = v7(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
grid%vv%v8(ips:ipe,jps:jpe,kps:kpe) = v8(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
grid%vv%v9(ips:ipe,jps:jpe,kps:kpe) = v9(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
grid%vv%v10(ips:ipe,jps:jpe,kps:kpe) = v10(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
end if
if ( use_cv_w ) then
grid%vv%v11(ips:ipe,jps:jpe,kps:kpe) = v11(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
end if
call da_vv_to_cv( grid%vv, grid%xp, be%cv_mz, be%ncv_mz, cv_size, cvt )
elseif ( use_inverse_squarerootb ) then
grid%vp%v1(ips:ipe,jps:jpe,kps:kpe) = v1(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
grid%vp%v2(ips:ipe,jps:jpe,kps:kpe) = v2(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
grid%vp%v3(ips:ipe,jps:jpe,kps:kpe) = v3(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
grid%vp%v4(ips:ipe,jps:jpe,kps:kpe) = v4(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
grid%vp%v5(ips:ipe,jps:jpe,kps:kpe) = v5(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
if ( cloud_cv_options >= 2 ) then
grid%vp%v6(ips:ipe,jps:jpe,kps:kpe) = v6(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
grid%vp%v7(ips:ipe,jps:jpe,kps:kpe) = v7(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
grid%vp%v8(ips:ipe,jps:jpe,kps:kpe) = v8(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
grid%vp%v9(ips:ipe,jps:jpe,kps:kpe) = v9(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
grid%vp%v10(ips:ipe,jps:jpe,kps:kpe) = v10(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
end if
if ( use_cv_w ) then ! vertical stagging +1?
grid%vp%v11(ips:ipe,jps:jpe,kps:kpe) = v11(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
end if
!call da_write_vp(grid,grid%vp,'vp_input.global ') ! to verify correctness
print '(/10X,"===> Use inverse transform of square-root B for outer-loop=",i2)', it
if ( cv_options == 3 ) then
Expand All @@ -639,6 +670,14 @@
deallocate( v3 )
deallocate( v4 )
deallocate( v5 )
if ( cloud_cv_options >= 2 ) then
deallocate( v6 )
deallocate( v7 )
deallocate( v8 )
deallocate( v9 )
deallocate( v10 )
end if
if ( use_cv_w ) deallocate( v11 )

close(vp_unit)
call da_free_unit(vp_unit)
Expand Down
6 changes: 3 additions & 3 deletions var/da/da_recursive_filter/da_transform_through_rf_inv.inc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,9 +33,9 @@ subroutine da_transform_through_rf_inv(grid, mz,rf_alpha, val, field, scaling)
if (trace_use_dull) call da_trace_entry("da_transform_through_rf_inv")

write (*,*) 'mz= ', mz
write (*,*) 'rf_alpha= ', rf_alpha
write (*,*) 'eigval= ', val
write (*,*) 'vert_corr=', vert_corr, ' vert_corr_1=', vert_corr_1
!write (*,*) 'rf_alpha= ', rf_alpha
!write (*,*) 'eigval= ', val
!write (*,*) 'vert_corr=', vert_corr, ' vert_corr_1=', vert_corr_1


rf_passes_over_two = rf_passes / 2
Expand Down
65 changes: 62 additions & 3 deletions var/da/da_setup_structures/da_write_vp.inc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ subroutine da_write_vp (grid,vp,filename)
! will be interpolated into higher resolution by offline program
! Method: based on da_write_increments.inc
! Author: Zhiquan (Jake) Liu, NCAR/MMM, 2015-09
! add cloud and w variables, 2017-07
!----------------------------------------------------------------------

implicit none
Expand All @@ -22,8 +23,9 @@ subroutine da_write_vp (grid,vp,filename)

!real, dimension(1:grid%xb%mix,1:grid%xb%mjy,1:grid%xb%mkz+1) :: wgbuf
real, dimension(:,:,:), allocatable :: v1_global, v2_global, &
v3_global, v4_global
real, dimension(:,:,:) , allocatable :: v5_global
v3_global, v4_global, v5_global
real, dimension(:,:,:), allocatable :: v6_global, v7_global, &
v8_global, v9_global, v10_global, v11_global
#endif

integer :: vp_unit, vp_local_unit
Expand All @@ -46,6 +48,16 @@ subroutine da_write_vp (grid,vp,filename)
allocate ( v3_global (1:ix,1:jy,1:kz))
allocate ( v4_global (1:ix,1:jy,1:kz))
allocate ( v5_global (1:ix,1:jy,1:kz))
if ( cloud_cv_options >= 2 ) then
allocate ( v6_global (1:ix,1:jy,1:kz))
allocate ( v7_global (1:ix,1:jy,1:kz))
allocate ( v8_global (1:ix,1:jy,1:kz))
allocate ( v9_global (1:ix,1:jy,1:kz))
allocate ( v10_global (1:ix,1:jy,1:kz))
end if
if ( use_cv_w ) then
allocate ( v11_global (1:ix,1:jy,1:kz))
end if

call da_patch_to_global(grid, vp % v1, gbuf) ! psi or u
if (rootproc) then
Expand All @@ -57,7 +69,6 @@ subroutine da_write_vp (grid,vp,filename)
v2_global(1:ix,1:jy,1:kz) = gbuf(1:ix,1:jy,1:kz)
end if

!call da_patch_to_global(grid, grid%xa % t, gbuf) ! t_u or t
call da_patch_to_global(grid, vp % v3, gbuf) ! t_u or t
if (rootproc) then
v3_global(1:ix,1:jy,1:kz) = gbuf(1:ix,1:jy,1:kz)
Expand All @@ -74,6 +85,40 @@ subroutine da_write_vp (grid,vp,filename)
v5_global(1:ix,1:jy,1:kz) = gbuf(1:ix,1:jy,1:kz)
end if

if ( cloud_cv_options >= 2 ) then
call da_patch_to_global(grid, vp % v6, gbuf) ! qcloud
if (rootproc) then
v6_global(1:ix,1:jy,1:kz) = gbuf(1:ix,1:jy,1:kz)
end if

call da_patch_to_global(grid, vp % v7, gbuf) ! qrain
if (rootproc) then
v7_global(1:ix,1:jy,1:kz) = gbuf(1:ix,1:jy,1:kz)
end if

call da_patch_to_global(grid, vp % v8, gbuf) ! qice
if (rootproc) then
v8_global(1:ix,1:jy,1:kz) = gbuf(1:ix,1:jy,1:kz)
end if

call da_patch_to_global(grid, vp % v9, gbuf) ! qsnow
if (rootproc) then
v9_global(1:ix,1:jy,1:kz) = gbuf(1:ix,1:jy,1:kz)
end if

call da_patch_to_global(grid, vp % v10, gbuf) ! qgraupel
if (rootproc) then
v10_global(1:ix,1:jy,1:kz) = gbuf(1:ix,1:jy,1:kz)
end if
end if ! cloud_cv_options

if ( use_cv_w ) then
call da_patch_to_global(grid, vp % v11, gbuf) ! w
if (rootproc) then
v11_global(1:ix,1:jy,1:kz) = gbuf(1:ix,1:jy,1:kz)
end if
end if

!write(unit=vpfile,fmt='(a,i4.4)') 'vp_',myproc
!call da_get_unit(vp_local_unit)
!open(unit=vp_local_unit, file=trim(vpfile), form='unformatted')
Expand Down Expand Up @@ -111,6 +156,12 @@ subroutine da_write_vp (grid,vp,filename)
write (unit=vp_unit) v1_global, v2_global, &
v3_global, v4_global, v5_global

if ( cloud_cv_options >= 2 ) then
write (unit=vp_unit) v6_global, v7_global, &
v8_global, v9_global, v10_global
end if
if ( use_cv_w ) write (unit=vp_unit) v11_global

close(vp_unit)
call da_free_unit(vp_unit)

Expand All @@ -122,6 +173,14 @@ subroutine da_write_vp (grid,vp,filename)
vp%v3(1:ix,1:jy,1:kz), &
vp%v4(1:ix,1:jy,1:kz), &
vp%v5(1:ix,1:jy,1)
if ( cloud_cv_options >= 2 ) then
write (unit=vp_unit) vp%v6(1:ix,1:jy,1:kz), &
vp%v7(1:ix,1:jy,1:kz), &
vp%v8(1:ix,1:jy,1:kz), &
vp%v9(1:ix,1:jy,1:kz), &
vp%v10(1:ix,1:jy,1:kz)
end if
if ( use_cv_w ) write (unit=vp_unit) vp%v11(1:ix,1:jy,1:kz)

close(vp_unit)
call da_free_unit(vp_unit)
Expand Down
182 changes: 4 additions & 178 deletions var/da/da_vtox_transforms/da_transform_vptox_inv.inc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,9 +24,10 @@ subroutine da_transform_vptox_inv(grid, vp, be, ep)
! [1] Add flow-dependent increments in control variable space (vp):
!---------------------------------------------------------------------------

if (be % ne > 0 .and. alphacv_method == alphacv_method_vp) then
call da_add_flow_dependence_vp(be % ne, ep, vp, its,ite, jts,jte, kts,kte)
end if
!if (be % ne > 0 .and. alphacv_method == alphacv_method_vp) then
! call da_add_flow_dependence_vp(be % ne, ep, vp, its,ite, jts,jte, kts,kte)
! call da_add_flow_dependence_vp_inv !!! ??
!end if

!--------------------------------------------------------------------------
! [2] Impose statistical balance constraints:
Expand Down Expand Up @@ -151,182 +152,7 @@ subroutine da_transform_vptox_inv(grid, vp, be, ep)
end if

end do
!$OMP END PARALLEL DO
!--------------------------------------------------------------------------
! [3] Transform to model variable space:
!--------------------------------------------------------------------------

!!#ifdef A2C
! if ((fg_format==fg_format_wrf_arw_regional .or. &
! fg_format==fg_format_wrf_arw_global ) .and. ide == ipe ) then
! ipe = ipe + 1
! ide = ide + 1
! end if
!
! if ((fg_format==fg_format_wrf_arw_regional .or. &
! fg_format==fg_format_wrf_arw_global ) .and. jde == jpe ) then
! jpe = jpe + 1
! jde = jde + 1
! end if
!!#endif

!!#ifdef DM_PARALLEL
!!#include "HALO_PSICHI_UV.inc"
!!#endif

!!#ifdef A2C
!! if ((fg_format==fg_format_wrf_arw_regional .or. &
! fg_format==fg_format_wrf_arw_global ) .and. ide == ipe ) then
! ipe = ipe - 1
! ide = ide - 1
! end if

! if ((fg_format==fg_format_wrf_arw_regional .or. &
! fg_format==fg_format_wrf_arw_global ) .and. jde == jpe ) then
! jpe = jpe - 1
! jde = jde - 1
! end if
!#endif

! Psi and chi to u and v:
! if ( cv_options == 5 .or. cv_options == 6 ) then
! call da_psichi_to_uv(vp % v1, vp % v2, grid%xb % coefx, &
! grid%xb % coefy , grid%xa % u, grid%xa % v)
! else if ( cv_options == 7 ) then
! grid%xa%u = vp%v1
! grid%xa%v = vp%v2
! end if

if ( (use_radarobs .and. use_radar_rf) .or. (use_rad .and. crtm_cloud).or. &
(use_radarobs .and. use_radar_rhv) .or. (use_radarobs .and. use_radar_rqv) .or. cloud_cv_options .ge. 2 .or. &
(grid%pseudo_var(1:1).eq.'q' .and. grid%pseudo_var(2:2).ne.' ') .or. &
(grid%pseudo_var(1:1).eq.'Q' .and. grid%pseudo_var(2:2).ne.' ') ) then

! if ( cloud_cv_options == 1 .and. use_3dvar_phy) then
! ! Pseudo RH --> Total water mixing ratio:
! !$OMP PARALLEL DO &
! !$OMP PRIVATE ( ij, i, j, k )
! do ij = 1 , grid%num_tiles
! do k = kts, kte
! do j = grid%j_start(ij), grid%j_end(ij)
! do i = its, ite
! grid%xa % qt(i,j,k) = vp%v4(i,j,k) * grid%xb%qs(i,j,k)
! enddo
! enddo
! enddo
! enddo
! !$OMP END PARALLEL DO
! end if
! if ( cloud_cv_options .ge. 2 ) then
! ! Pseudo RH --> Water vapor mixing ratio:
! !$OMP PARALLEL DO &
! !$OMP PRIVATE ( ij, i, j, k )
! do ij = 1 , grid%num_tiles
! do k = kts, kte
! do j = grid%j_start(ij), grid%j_end(ij)
! do i = its, ite
! grid%xa % q(i,j,k) = vp%v4(i,j,k) * grid%xb%qs(i,j,k)
! enddo
! enddo
! enddo
! enddo
! !$OMP END PARALLEL DO
#ifdef CLOUD_CV
!qcloud
!$OMP PARALLEL DO &
!$OMP PRIVATE ( ij, i, j, k )
do ij = 1 , grid%num_tiles
do k = kts, kte
do j = grid%j_start(ij), grid%j_end(ij)
do i = its, ite
vp%v6(i,j,k) = grid%xa % qcw(i,j,k)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
!qrain
!$OMP PARALLEL DO &
!$OMP PRIVATE ( ij, i, j, k )
do ij = 1 , grid%num_tiles
do k = kts, kte
do j = grid%j_start(ij), grid%j_end(ij)
do i = its, ite
vp%v7(i,j,k) = grid%xa % qrn(i,j,k)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
!qice
!$OMP PARALLEL DO &
!$OMP PRIVATE ( ij, i, j, k )
do ij = 1 , grid%num_tiles
do k = kts, kte
do j = grid%j_start(ij), grid%j_end(ij)
do i = its, ite
vp%v8(i,j,k) = grid%xa % qci(i,j,k)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
!qsnow
!$OMP PARALLEL DO &
!$OMP PRIVATE ( ij, i, j, k )
do ij = 1 , grid%num_tiles
do k = kts, kte
do j = grid%j_start(ij), grid%j_end(ij)
do i = its, ite
vp%v9(i,j,k) = grid%xa % qsn(i,j,k)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
!qgraupel
!$OMP PARALLEL DO &
!$OMP PRIVATE ( ij, i, j, k )
do ij = 1 , grid%num_tiles
do k = kts, kte
do j = grid%j_start(ij), grid%j_end(ij)
do i = its, ite
vp%v10(i,j,k) = grid%xa % qgr(i,j,k)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
!vertical velocity
!$OMP PARALLEL DO &
!$OMP PRIVATE ( ij, i, j, k )
do ij = 1 , grid%num_tiles
do k = kts, kte
do j = grid%j_start(ij), grid%j_end(ij)
do i = its, ite
vp%v11(i,j,k) = grid%xa % w(i,j,k)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
#endif
! end if
!else ! no rf or cloud radiance
! ! Pseudo RH --> Water vapor mixing ratio:
! !$OMP PARALLEL DO &
! !$OMP PRIVATE ( ij, i, j, k )
! do ij = 1 , grid%num_tiles
! do k = kts, kte
! do j = grid%j_start(ij), grid%j_end(ij)
! do i = its, ite
! grid%xa % q(i,j,k) = vp%v4(i,j,k) * grid%xb%qs(i,j,k)
! enddo
! enddo
! enddo
! enddo
! !$OMP END PARALLEL DO
end if ! RF or Radiance
!---------------------------------------------------------------------------
! [4] Add flow-dependent increments in model space (grid%xa):
!---------------------------------------------------------------------------
Expand Down
Loading

0 comments on commit 606ac0e

Please sign in to comment.