wrfv2_fire/var/da/da_minimisation/da_calculate_j.inc

   1 subroutine da_calculate_j(it, iter, cv_size, cv_size_jb, cv_size_je, cv_size_jp, &
   2                            xbx, be, iv, xhat, cv, &
   3                            re, y, j, grid, config_flags                     )
   4
   5    !---------------------------------------------------------------------------
   6    ! Purpose: Initialises the Y-array
   7    !---------------------------------------------------------------------------
   8
   9    implicit none
  10
  11    integer, intent(in)                :: it     ! external iteration #.
  12    integer, intent(in)                :: iter   ! internal iteration #.
  13    integer, intent(in)                :: cv_size    ! Total cv size.
  14    integer, intent(in)                :: cv_size_jb ! Jb cv size.
  15    integer, intent(in)                :: cv_size_je ! Je cv size.
  16    integer, intent(in)                :: cv_size_jp ! Jp cv size.
  17    type (xbx_type),intent(in)         :: xbx    ! For header & non-grid arrays.
  18    type (be_type), intent(in)         :: be     ! background error structure.
  19    type (iv_type), intent(inout)      :: iv     ! innovation vector (o-b).
  20    real, intent(in)                   :: xhat(1:cv_size) ! control variables.
  21    real, intent(in)                   :: cv(1:cv_size)   ! control variables.
  22    type (y_type) , intent(inout)      :: re     ! residual vector (o-a).
  23    type (y_type) , intent(inout)      :: y      ! y = H(x_inc).
  24    type (j_type) , intent(out)        :: j      ! cost function j
  25
  26    type(domain), intent(inout)  :: grid
  27    type(grid_config_rec_type), intent(inout) :: config_flags
  28
  29    integer          :: je_start, je_end             ! Start/end indices of Je.
  30    real             :: jo_partial                   ! jo for this processor
  31    type (y_type)    :: jo_grad_y ! Grad_y(jo)
  32    real             :: cv_xhat_jb(cv_size_jb), cv_xhat_je(cv_size_je)
  33
  34    integer          :: ndynopt, is, ie, js, je, ks, ke
  35    real             :: dtemp1x
  36
  37    ! Variables for VarBC background constraint
  38    real :: cv_xhat_jp(cv_size_jp), xhat_jp(cv_size_jp)    ! Jp control variable.
  39    integer                           :: jp_start, jp_end       ! Start/end indices of Jp.
  40    integer                           :: inst, ichan, npred, ipred, id
  41    real                              :: bgerr, gnorm_jp
  42
  43    integer                           :: n, cldtoplevel(1), icld, nclouds, ncv, minlev_cld
  44    real                              :: js_local
  45    real, allocatable                 :: cc(:)
  46
  47    if (trace_use) call da_trace_entry("da_calculate_j")
  48
  49    je_start = cv_size_jb + 1
  50    je_end   = cv_size_jb + cv_size_je
  51
  52    jp_start = cv_size_jb + cv_size_je + 1
  53    jp_end   = cv_size_jb + cv_size_je + cv_size_jp
  54
  55    call da_allocate_y(iv, jo_grad_y)
  56
  57    !-------------------------------------------------------------------------
  58    ! [1.0] calculate jo:
  59    !-------------------------------------------------------------------------
  60
  61    ! [1.1] transform from control variable to model grid space:
  62
  63    if (iter > 0) &
  64       call da_transform_vtoy(cv_size, be, grid%ep, xhat, iv, grid%vp, grid%vv,&
  65                               xbx, y, &
  66                               grid, config_flags                      )
  67
  68    ! [1.2] compute residual (o-a) = (o-b) - h x~
  69
  70    call da_calculate_residual(iv, y, re)
  71
  72    ! [1.3] calculate jo:
  73
  74    call da_jo_and_grady(iv, re, jo_partial, j % jo, jo_grad_y)
  75
  76    if (test_dm_exact) then
  77       ! jo_partial has been already summed at lower level
  78       j % jo % total = jo_partial
  79    else
  80       j % jo % total = wrf_dm_sum_real(jo_partial)
  81    end if
  82
  83    ! [1.4] calculate jc-dfi:
  84
  85    j % jc = 0.0
  86
  87    if ( var4d .and. grid%jcdfi_io .and. iter > 0 ) then
  88
  89       ndynopt      = grid%dyn_opt
  90       grid%dyn_opt = DYN_EM_TL
  91       call nl_set_dyn_opt (1 , DYN_EM_TL)
  92
  93       call da_med_initialdata_input(grid , config_flags, 'tldf')
  94
  95       grid%dyn_opt = ndynopt
  96       call nl_set_dyn_opt (1 , DYN_EM)
  97
  98       is = grid%xp%ips
  99       ie = grid%xp%ipe
 100       if ( ie == grid%xp%ide ) ie = ie + 1
 101       js = grid%xp%jps
 102       je = grid%xp%jpe
 103       if ( je == grid%xp%jde ) je = je + 1
 104       ks = grid%xp%kps
 105       ke = grid%xp%kpe
 106
 107       j % jc = j % jc + 0.5*grid%jcdfi_gama*da_dot( (ie-is+1)*(je-js+1)*(ke-ks+1), grid%g_u_2(is:ie,js:je,ks:ke), grid%g_u_2(is:ie,js:je,ks:ke) )
 108       j % jc = j % jc + 0.5*grid%jcdfi_gama*da_dot( (ie-is+1)*(je-js+1)*(ke-ks+1), grid%g_v_2(is:ie,js:je,ks:ke), grid%g_v_2(is:ie,js:je,ks:ke) )
 109       j % jc = j % jc + 0.5*grid%jcdfi_gama*da_dot( (ie-is+1)*(je-js+1)*(ke-ks+1), grid%g_t_2(is:ie,js:je,ks:ke), grid%g_t_2(is:ie,js:je,ks:ke) )
 110       j % jc = j % jc + 0.5*grid%jcdfi_gama*da_dot( (ie-is+1)*(je-js+1), grid%g_mu_2(is:ie,js:je), grid%g_mu_2(is:ie,js:je) )
 111       j % jc = j % jc + 0.5*grid%jcdfi_gama*da_dot( (ie-is+1)*(je-js+1)*(ke-ks+1), grid%g_moist(is:ie,js:je,ks:ke,P_G_QV), grid%g_moist(is:ie,js:je,ks:ke,P_G_QV) )
 112
 113       dtemp1x = j % jc
 114       ! summation across processors:
 115       j % jc  = wrf_dm_sum_real(dtemp1x)
 116
 117    end if
 118
 119    !-------------------------------------------------------------------------
 120    ! [2.0] calculate jb:
 121    !-------------------------------------------------------------------------
 122
 123    cv_xhat_jb(1:cv_size_jb) = cv(1:cv_size_jb) + xhat(1:cv_size_jb)
 124
 125    j % jb = 0.5 * da_dot_cv(cv_size_jb, cv_size_domain_jb, &
 126        cv_xhat_jb, cv_xhat_jb, grid, &
 127        (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
 128
 129    j % jb = jb_factor * j % jb
 130
 131    !-------------------------------------------------------------------------
 132    ! [3.0] calculate je:
 133    !-------------------------------------------------------------------------
 134
 135    j % je = 0.0
 136    if (be % ne > 0) then
 137       cv_xhat_je(1:cv_size_je) = cv(je_start:je_end) + xhat(je_start:je_end)
 138
 139       j % je = 0.5 * da_dot_cv(cv_size_je, cv_size_domain_je, &
 140           cv_xhat_je, cv_xhat_je, grid, &
 141           (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
 142    end if
 143    j % je = je_factor * j % je
 144
 145    !-------------------------------------------------------------------------
 146    ! [4.0] calculate jp:
 147    !-------------------------------------------------------------------------
 148    j % jp = 0.0
 149 #if defined(RTTOV) || defined(CRTM)
 150    if (use_varbc) then
 151       do inst = 1, iv % num_inst
 152          do ichan = 1, iv%instid(inst)%nchan
 153             if (satinfo(inst)%iuse(ichan) <= 0) cycle
 154             npred    = iv%instid(inst)%varbc(ichan)%npred
 155             if (npred <= 0) cycle               !! VarBC channels only
 156             do ipred = 1, npred
 157                id     = iv%instid(inst)%varbc(ichan)%index(ipred)
 158                bgerr  = iv%instid(inst)%varbc(ichan)%bgerr(ipred)
 159                if (bgerr > 0.0) &
 160                   cv_xhat_jp(id-jp_start+1) = (1/sqrt(bgerr)) * &
 161                      SUM( (cv(id)+xhat(id)) * &
 162                      iv%instid(inst)%varbc(ichan)%vtox(ipred,1:npred))
 163             end do
 164          end do
 165       end do
 166       j % jp = 0.5 * da_dot(cv_size_jp, cv_xhat_jp, cv_xhat_jp)
 167    end if
 168 #endif
 169
 170    !-------------------------------------------------------------------------
 171    ! [5.0] calculate js:
 172    !-------------------------------------------------------------------------
 173    j % js = 0.0
 174    if (ANY(use_satcv)) then
 175       do inst = 1, iv % num_inst
 176          do n = iv%instid(inst)%info%n1, iv%instid(inst)%info%n2 ! loop for pixel
 177          ! Skin Temperature
 178          !-----------------
 179             if (use_satcv(1)) then
 180 !               j % js = j % js + 0.5 * xhat(iv%instid(inst)%cv_index(n)%ts) **2
 181
 182 !              !!! Super-TMP dump of Tskin increment for plotting purposes
 183 !               if (iter > 0) iv%instid(inst)%tb_xb(1,n)  = xhat(iv%instid(inst)%cv_index(n)%ts)
 184             end if
 185
 186          ! Cloud cover(s)
 187          !---------------
 188             if (use_satcv(2)) then
 189             j % js = j % js + 0.5 * SUM( xhat(iv%instid(inst)%cv_index(n)%cc) **2)
 190
 191             j % js = j % js + 0.5 * SUM( (10.0 * xhat(iv%instid(inst)%cv_index(n)%cc)) **2,      &
 192                                           MASK = xhat(iv%instid(inst)%cv_index(n)%cc) < 0.0 .or. &
 193                                                  xhat(iv%instid(inst)%cv_index(n)%cc) > 1.0 )
 194
 195                if (iter > 0) then
 196                   nclouds = iv%instid(inst)%cv_index(n)%nclouds
 197                   ncv     = iv%instid(inst)%cv_index(n)%ncv
 198                   allocate(cc(nclouds))
 199
 200                   cc = xhat(iv%instid(inst)%cv_index(n)%cc)
 201                !---------------------------------------------------------------
 202                ! Change of variable (preconditioning)
 203                !---------------------------------------------------------------
 204 !                 do icld = 1, nclouds
 205 !                    cc(icld) = SUM( xhat(iv%instid(inst)%cv_index(n)%cc) * &
 206 !                                       iv%instid(inst)%cv_index(n)%vtox(icld,1:ncv) )
 207 !                 end do
 208
 209                   if (use_satcv(1)) then
 210                      write (*, '(i6,100F8.2)')n,xhat(iv%instid(inst)%cv_index(n)%ts), SUM(cc)*100, cc*100
 211                   else
 212                      write (*, '(i6,100F8.2)')n,SUM(cc)*100, cc*100
 213                   end if
 214
 215 !                 !!! Super-TMP dump of Cloud Cover increment for plotting purposes
 216 !                  iv%instid(inst)%tb_inv(1,n) = SUM(cc)*100.0
 217 !
 218 !                 !!! Super-TMP dump of Cloud Top Pressure for plotting purposes
 219 !                 minlev_cld = 5
 220 !                 if (ANY(cc(minlev_cld:nclouds) > 0.01)) then
 221 !                    cldtoplevel = MINLOC(cc(minlev_cld:nclouds), MASK = cc(minlev_cld:nclouds) > 0.01)
 222 !                 else
 223 !                    cldtoplevel = nclouds
 224 !                 end if
 225 !                 cldtoplevel = cldtoplevel + kte - nclouds !!!+ minlev_cld
 226 !!                  if (rtm_option == rtm_option_rttov) then
 227 !!                     re%instid(inst)%tb(1,n) = coefs(inst)%ref_prfl_p(cldtoplevel(1))
 228 !!                  elseif (rtm_option == rtm_option_crtm) then
 229 !                     re%instid(inst)%tb(1,n) = iv%instid(inst)%pm(cldtoplevel(1),n)
 230 !!                  end if
 231
 232                   deallocate(cc)
 233                end if
 234             end if
 235          end do
 236       end do
 237       js_local = j % js
 238       ! summation across processors:
 239       j % js = wrf_dm_sum_real(js_local)
 240    end if
 241
 242    !-------------------------------------------------------------------------
 243    ! [6.0] calculate total cost function j = jo + jb + jc + je + jp + js:
 244    !-------------------------------------------------------------------------
 245
 246    j % total = j % jb + j % jo % total + j % je + j % jp + j % js
 247    if (grid%jcdfi_use) j % total = j % total  + j % jc
 248
 249    !-------------------------------------------------------------------------
 250    ! [7.0] write cost function:
 251    !-------------------------------------------------------------------------
 252    if (it == 1 .and. iter == 0 .and. rootproc) then
 253       write(unit=cost_unit,fmt='(a)')'Outer    EPS     Inner      J           Jb       Jo           Jc         Je         Jp         Js'
 254       write(unit=cost_unit,fmt='(a)')'Iter             Iter                            '
 255    end if
 256
 257    if (rootproc) then
 258       write(unit=cost_unit,fmt='(2x,i2,1x,e10.3,2x,i4,7(1x,f10.3))') &
 259          it, EPS(it), iter, j % total, j % jb, j % jo % total, j % jc, j % je, j % jp, j%js
 260    end if
 261
 262    if (trace_use) call da_trace_exit("da_calculate_j")
 263
 264 end subroutine da_calculate_j
 265