Add the new capability for MRI-4DVar

Registry/registry.var

@@ -264,9 +264,10 @@ rconfig   integer   report_end              namelist,wrfvar5  1  10000000 - "rep
 rconfig   integer   tovs_start              namelist,wrfvar5  1  1        - "tovs_start"               ""  ""
 rconfig   integer   tovs_end                namelist,wrfvar5  1  10000000 - "tovs_end"                 ""  ""
 rconfig   logical   gpsref_thinning         namelist,wrfvar5  1  .false.  - "gpsref_thinning"          ""  ""
-rconfig   logical   outer_loop_restart      namelist,wrfvar6  1  .false.  - "outer_loop_restart"       ""  ""
 rconfig   integer   max_ext_its             namelist,wrfvar6  1  1        - "max_ext_its"              ""  ""
 rconfig   integer   ntmax                   namelist,wrfvar6  max_outer_iterations  75      - "ntmax"                    ""  ""
+rconfig   logical   use_inverse_squarerootb namelist,wrfvar6  1  .false.  - "use_inverse_squarerootb"       ""  ""
+rconfig   logical   use_interpolate_cvt     namelist,wrfvar6  1  .false.  - "use_interpolate_cvt"       ""  ""
 rconfig   integer   nsave                   namelist,wrfvar6  1  4        - "nsave"                    ""  ""
 rconfig   integer   write_interval          namelist,wrfvar6  1  5        - "write_interval"           ""  ""
 rconfig   real      eps                     namelist,wrfvar6  max_outer_iterations 0.01     - "eps"                      ""  ""

var/da/da_main/da_solve.inc

@@ -45,6 +45,8 @@
    type(x_type) :: shuffle
 
    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
@@ -53,8 +55,9 @@
    real,    external :: nest_loc_of_cg  ! from share/interp_fcn.F
    integer, external :: compute_CGLL    ! from share/interp_fcn.F
 
-   integer   :: cvt_unit, iost
-   character(len=8) :: cvtfile
+   integer   :: vp_unit, iost
+   character(len=13) :: vpfile ! vp_input.0001
+   integer   :: i1,i2,i3,i4,i5,i6
    logical :: ex
 
    character(len=10) :: this_time
@@ -630,42 +633,135 @@
          ! Done with randomcv.
          ! Set the following to skip some code to get to the deallocation part.
          max_ext_its = 0
-         outer_loop_restart = .false.
       end if  !anal_type_randomcv
 
-      if ( outer_loop_restart ) then
-         !call da_get_unit(cvt_unit)
-         cvt_unit=600
+! mri-4dvar: if multi_inc /= 2: run normal 3D/4D-Var
+!------------------------------------------------------------------------
+    ! cvt is outer loop control variable, it is zero for the first outer loop,
+    ! but non-zero from the second outer loop in normal 3d/4dvar.
+    ! for MRI-4DVar, vp from the previous outer loop needs to be read in, 
+    ! then perform the inverse transform to derive cvt
+    !-----------------------------------------------------
+      call da_initialize_cv (cv_size, cvt)
+      call da_zero_vp_type (grid%vp)
+      call da_zero_vp_type (grid%vv)
+
+      if ( multi_inc == 2 ) then
          if ( max_ext_its > 1 ) then
            max_ext_its=1
-           write(unit=message(1),fmt='(a)') "Re-set max_ext_its = 1 for outer_loop_restart"
+           write(unit=message(1),fmt='(a)') "Re-set max_ext_its = 1 for multi_inc==2"
            call da_message(message(1:1))
          end if
-         write(unit=cvtfile,fmt='(a,i4.4)') 'cvt_',myproc
-         inquire(file=trim(cvtfile), exist=ex)
+
+      ! read vp files for different PEs
+      !----------------------------------
+         write(unit=vpfile,fmt='(a,i4.4)') 'vp_input.',myproc
+         inquire(file=trim(vpfile), exist=ex)
          if ( ex ) then
-           open(unit=cvt_unit,file=trim(cvtfile),iostat=iost,form='UNFORMATTED',status='OLD')
+           call da_get_unit(vp_unit)
+           open(unit=vp_unit,file=trim(vpfile),iostat=iost,form='UNFORMATTED',status='OLD')
            if (iost /= 0) then
              write(unit=message(1),fmt='(A,I5,A)') &
-                "Error ",iost," opening cvt file "//trim(cvtfile)
+                "Error ",iost," opening vp file "//trim(vpfile)
              call da_error(__FILE__,__LINE__,message(1:1))
            end if
-           write(unit=message(1),fmt='(a)') 'Reading cvt from : '//trim(cvtfile) 
+           if ( use_interpolate_cvt ) then         ! works for CV3?, 3D RF
+               write(unit=message(1),fmt='(a)') 'Reading vv from : '//trim(vpfile) 
+           elseif ( use_inverse_squarerootb ) then ! works for CV5,6,7, vertical EOF
+               write(unit=message(1),fmt='(a)') 'Reading vp from : '//trim(vpfile)
+           end if 
            call da_message(message(1:1))
-           read(cvt_unit) cvt
-           close(cvt_unit)
+           read(vp_unit) i1, i2, i3, i4, i5, i6 ! read dimension of patch for current processor
+           allocate( v1(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1) )
+           allocate( v2(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1) )
+           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
+                write(unit=message(1),fmt='(A,I5,A)') &
+                   "Error: inverse U transform not for cv_options = 3"
+                call da_error(__FILE__,__LINE__,message(1:1))
+             end if
+             call da_transform_vtox_inv (grid,be%cv%size_jb,xbx,be,grid%ep,cvt(1:be%cv%size_jb),grid%vv,grid%vp)
+           end if
+
+           deallocate( v1 )
+           deallocate( v2 )
+           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)
+
          else
-           write(unit=message(1),fmt='(a)') "cvt file '"//trim(cvtfile)//"' does not exists, initializing cvt."
+           write(unit=message(1),fmt='(a)') "vp files '"//trim(vpfile)//"' does not exists, initiallizing cvt."
            call da_message(message(1:1))
-           call da_initialize_cv (cv_size, cvt)
+           call da_initialize_cv (cv_size, cvt) ! perhaps better use da_error
          end if
-      else
-         call da_initialize_cv (cv_size, cvt)
+
+         call da_zero_vp_type (grid%vv)
+         call da_zero_vp_type (grid%vp)
+
       end if
+! mri-4dvar -------------------------------------------
 
-      call da_zero_vp_type (grid%vv)
-      call da_zero_vp_type (grid%vp)
-
       if ( var4d ) then
 #ifdef VAR4D
          call da_zero_vp_type (grid%vv6)
@@ -698,6 +794,24 @@
 
       call da_initialize_cv (cv_size, xhat)
 
+! mri-4dvar----------------------
+! Apply inverse transform of squareroot(B) for full-resolution non-stop Var
+! from 2nd outer loop, this is to check correctness of inverse U transform
+! does not apply this setting for real world application
+!-----------------------------
+      if (multi_inc == 0 .and. it > 1 .and. use_inverse_squarerootb .and. cv_options /= 3) then
+         print '(/10X,"===> Use inverse transform of square-root B for outer-loop=",i2)', it
+         call da_transform_vtox_inv (grid,be%cv%size_jb,xbx,be,grid%ep,cvt(1:be%cv%size_jb),grid%vv,grid%vp)
+      endif
+
+! Reinitialize cvt=0 for full-resolution non-stop Var for each loop
+!------another option not tested --------------
+      if (multi_inc == 0 .and. it > 1 .and. use_interpolate_cvt) then
+         print '(/10X,"===> Reinitialize cvt as zeros for outer loop ",i2)', it
+         call da_initialize_cv (cv_size, cvt)
+      endif
+! mri-4dvar------------------------
+
       ! [8.1] Calculate nonlinear model trajectory 
 
 !     if (var4d .and. multi_inc /= 2 ) then
@@ -869,18 +983,9 @@
 
       ! Update outer-loop control variable
       cvt = cvt + xhat
-      if ( outer_loop_restart ) then
-        open(unit=cvt_unit,status='unknown',file=trim(cvtfile),iostat=iost,form='UNFORMATTED')
-        if (iost /= 0) then
-           write(unit=message(1),fmt='(A,I5,A)') &
-             "Error ",iost," opening cvt file "//trim(cvtfile)
-          call da_error(__FILE__,__LINE__,message(1:1))
-        end if
-        write(unit=message(1),fmt='(a)') 'Writing cvt to : '//trim(cvtfile) 
-        call da_message(message(1:1))
-        write(cvt_unit) cvt
-        close(cvt_unit)
-        !call da_free_unit(cvt_unit)
+      if ( multi_inc == 2 .and. use_interpolate_cvt ) then ! obsolete option
+        call da_cv_to_vv( cv_size, cvt, be%cv_mz, be%ncv_mz, grid%vv )
+        call da_write_vp(grid,grid%vv,'vp_output.global') ! wrtie vv to vp file
       end if
       !------------------------------------------------------------------------
 
@@ -905,6 +1010,13 @@
          call da_transform_vtox (grid,be%cv%size_jb,xbx,be,grid%ep,xhat(1:be%cv%size_jb),grid%vv,grid%vp)
          call da_transform_vpatox (grid,be,grid%ep,grid%vp)
       endif
+
+! mri-4dvar--------------------------
+      if (multi_inc == 2 .and. use_inverse_squarerootb) then
+         call da_write_vp(grid,grid%vp,'vp_output.global') ! write vp to vp file
+      end if
+! mri-4dvar--------------------------
+
       call da_transform_xtoxa (grid)
 
       ! [8.6] Only when use_radarobs = .false. and calc_w_increment =.true.,
@@ -1105,11 +1217,5 @@
 
    if (trace_use) call da_trace_exit ("da_solve")
 
-
-contains
-
-#include "da_solve_init.inc"
-#include "da_solve_dual_res_init.inc"
-
 end subroutine da_solve
 

var/da/da_main/da_wrfvar_top.f90

@@ -55,7 +55,8 @@ module da_wrfvar_top
    use da_obs, only : da_transform_xtoy_adj 
    use da_obs_io, only : da_write_filtered_obs, da_write_obs, da_final_write_obs , &
                          da_write_obs_etkf, da_write_modified_filtered_obs
-   use da_par_util, only : da_system,da_copy_tile_dims,da_copy_dims
+   use da_par_util, only : da_system,da_copy_tile_dims,da_copy_dims, &
+                           da_vv_to_cv, da_cv_to_vv
    use da_physics, only : da_uvprho_to_w_lin
 #if defined (CRTM) || defined (RTTOV)
    use da_radiance, only : da_deallocate_radiance
@@ -65,7 +66,7 @@ module da_wrfvar_top
    use da_varbc, only : da_varbc_init,da_varbc_update
 #endif
    use da_reporting, only : message, da_warning, da_error, da_message
-   use da_setup_structures, only : da_setup_obs_structures, &
+   use da_setup_structures, only : da_setup_obs_structures, da_write_vp, &
       da_setup_background_errors,da_setup_flow_predictors, &
       da_setup_cv, da_scale_background_errors, da_scale_background_errors_cv3
    use da_setup_structures, only : da_setup_flow_predictors_para_read_opt1
@@ -76,7 +77,8 @@ module da_wrfvar_top
    use da_transfer_model, only : da_transfer_xatoanalysis,da_setup_firstguess, &
        da_transfer_wrftltoxa_adj
    use da_vtox_transforms, only : da_transform_vtox, da_transform_xtoxa, &
-      da_transform_xtoxa_adj, da_copy_xa, da_add_xa, da_transform_vpatox
+      da_transform_xtoxa_adj, da_copy_xa, da_add_xa, da_transform_vpatox, &
+      da_transform_vtox_inv
    use da_wrfvar_io, only : da_med_initialdata_input, da_update_firstguess
    use da_tools, only : da_set_randomcv, da_get_julian_time
 
@@ -150,5 +152,7 @@ module da_wrfvar_top
 #include "da_wrfvar_interface.inc"
 #include "da_wrfvar_finalize.inc"
 #include "da_solve.inc"
+#include "da_solve_init.inc"
+#include "da_solve_dual_res_init.inc"
 
 end module da_wrfvar_top

var/da/da_minimisation/da_get_innov_vector.inc

@@ -2,9 +2,9 @@ subroutine da_get_innov_vector (it, num_qcstat_conv, ob, iv, grid, config_flags)
 
 
    !-----------------------------------------------------------------------
-   ! Purpose: TBD
-   !    Updated for Analysis on Arakawa-C grid
-   !    Author: Syed RH Rizvi,  MMM/ESSL/NCAR,  Date: 10/22/2008
+   ! Purpose: driver routine for getting innovation vectors
+   ! History:$
+   !    10/22/2008 - Updated for Analysis on Arakawa-C grid (Syed RH Rizvi, NCAR)
    !-----------------------------------------------------------------------
 
    implicit none
@@ -145,7 +145,7 @@ subroutine da_get_innov_vector (it, num_qcstat_conv, ob, iv, grid, config_flags)
          call da_get_innov_vector_ssmt2     (it, num_qcstat_conv,grid, ob, iv)
       if (iv%info(satem)%nlocal     > 0) &
          call da_get_innov_vector_satem     (it, num_qcstat_conv,grid, ob, iv)
-      if (iv%info(radar)%nlocal     > 0) &
+      if (iv%info(radar)%nlocal     >= 0 .and. use_radarobs) &
          call da_get_innov_vector_radar     (it, grid, ob, iv)
       if (iv%info(qscat)%nlocal     > 0) &
          call da_get_innov_vector_qscat     (it, num_qcstat_conv,grid, ob, iv)
@@ -168,7 +168,7 @@ subroutine da_get_innov_vector (it, num_qcstat_conv, ob, iv, grid, config_flags)
       end if
       if (iv%info(airsr)%nlocal     > 0) &
       call da_get_innov_vector_airsr    (it,num_qcstat_conv, grid, ob, iv)
-
+ 
    !----------------------------------------------
    ! [5]  write out iv in ascii format
    !-----------------------------------------------
@@ -255,6 +255,18 @@ subroutine da_get_innov_vector (it, num_qcstat_conv, ob, iv, grid, config_flags)
       call domain_clock_set (grid, time_step_seconds=time_step_seconds)
       call domain_clockprint(150, grid, 'get CurrTime from clock,')
    end if
+#if defined(CRTM) || defined(RTTOV)
+   !----------------------------------------------
+   ! write out or read in radiance iv for multi in binary format
+   !-----------------------------------------------
+   if (use_rad) then
+      if ( multi_inc == 1 ) then
+         call da_write_iv_rad_for_multi_inc(it,ob,iv)
+      elseif ( multi_inc == 2 ) then
+         call da_read_iv_rad_for_multi_inc(it,ob,iv)
+      end if
+   end if
+#endif
 
    if ( multi_inc == 1 ) then
 #ifdef DM_PARALLEL

var/da/da_minimisation/da_minimisation.f90

@@ -54,7 +54,7 @@ module da_minimisation
       use_satcv, sensitivity_option, print_detail_outerloop, adj_sens, filename_len, &
       ims, ime, jms, jme, kms, kme, ips, ipe, jps, jpe, kps, kpe, fgat_rain_flags, var4d_bin_rain, freeze_varbc, &
       use_wpec, wpec_factor, use_4denvar, anal_type_hybrid_dual_res, alphacv_method, alphacv_method_xa, &
-      write_detail_grad_fn, pseudo_uvtpq, lanczos_ep_filename, use_divc, divc_factor, &
+      write_detail_grad_fn, pseudo_uvtpq, lanczos_ep_filename, use_divc, divc_factor, use_radarobs, &
       cloud_cv_options, use_cv_w, var_scaling6, var_scaling7, var_scaling8, var_scaling9, &
       var_scaling10, var_scaling11, &
       write_gts_omb_oma, write_unpert_obs, write_rej_obs_conv, pseudo_time, &
@@ -121,6 +121,7 @@ module da_minimisation
       da_get_innov_vector_radiance, satinfo
    use da_radiance1, only : da_ao_stats_rad,da_oi_stats_rad, &
       da_write_iv_rad_ascii,da_residual_rad,da_jo_and_grady_rad, &
+      da_write_iv_rad_for_multi_inc,da_read_iv_rad_for_multi_inc, &
       da_biasprep, da_qc_rad
 #endif
    use da_radar, only :  da_calculate_grady_radar, da_ao_stats_radar, &

var/da/da_obs_io/da_obs_io.f90

@@ -32,7 +32,7 @@ module da_obs_io
       thin_conv, thin_conv_ascii, lsac_nh_step, lsac_nv_step, lsac_nv_start, lsac_print_details, &
       lsac_use_u, lsac_use_v, lsac_use_t, lsac_use_q, lsac_u_error, lsac_v_error, lsac_t_error, lsac_q_error, &
       gpsro_drift, max_gpseph_input, use_gpsephobs, gpseph, gpseph_loadbalance, kds, kde, kts, kte, &
-      use_radar_rhv, use_radar_rqv
+      use_radar_rhv, use_radar_rqv, use_radar_rf, use_radar_rv, multi_inc
 
    use da_wrf_interfaces, only : wrf_dm_bcast_integer, wrf_dm_bcast_real
    use da_define_structures, only : iv_type, multi_level_type, multi_level_type_BUFR, &