subroutine da_calculate_gradj(it, iter, cv_size, cv_size_jb, cv_size_je, cv_size_jp, & cv_size_jl, cv_size_jt, xbx, be, iv, cv, y, grad, grid, config_flags, re ) !--------------------------------------------------------------------------- ! Purpose: Calculates the gradient of the cost function w/r to cv ! ! Called from da_minimise_cg (or da_minimise_lz) ! ! History: 12/12/08 - Creation from da_calculate_j (Tom Auligne) ! !--------------------------------------------------------------------------- implicit none integer, intent(in) :: it ! external iteration #. integer, intent(in) :: iter ! internal iteration #. integer, intent(in) :: cv_size ! Total cv size. integer, intent(in) :: cv_size_jb, cv_size_je, cv_size_jp, cv_size_jl, cv_size_jt type (xbx_type),intent(inout) :: xbx ! For header & non-grid arrays. type (be_type), intent(in) :: be ! background error structure. type (iv_type), intent(inout) :: iv ! innovation vector (o-b). real, intent(in) :: cv (1:cv_size) ! control variables. type (y_type), intent(inout) :: y real, intent(out) :: grad(cv_size) ! gradient of cost function type (y_type), optional, intent(inout) :: re ! residual vector (o-a). type(domain), intent(inout) :: grid type(grid_config_rec_type), intent(inout) :: config_flags integer :: je_start, je_end ! Start/end indices of Je. integer :: jl_start, jl_end ! Start/end indices of Jl. integer :: jt_start, jt_end, iphase real :: jo_partial ! jo for this processor type (y_type) :: jo_grad_y ! Grad_y(jo) real :: grad_jo(cv_size) real :: grad_jb(cv_size) real :: grad_je(cv_size) real :: grad_jd(cv_size) real :: grad_jp(cv_size) real :: grad_js(cv_size) real :: grad_jl(cv_size) real :: grad_jm(cv_size) real :: grad_jt(cv_size) real :: gnorm_j, gnorm_jo, gnorm_jb, gnorm_je, gnorm_jd, gnorm_jp, gnorm_js, gnorm_jl, gnorm_jt real :: gnorm_jm logical :: jcdf_flag real :: inc_div(ims:ime,jms:jme,kms:kme) ! Temp storage ! Variables for VarBC background constraint integer :: jp_start, jp_end ! Start/end indices of Jp. integer :: inst, ichan, npred, ipred, id real :: bgerr integer :: n if (trace_use) call da_trace_entry("da_calculate_gradj") !------------------------------------------------------------------------- ! [0.0] initialization: !------------------------------------------------------------------------- je_start = cv_size_jb + 1 je_end = cv_size_jb + cv_size_je jp_start = cv_size_jb + cv_size_je + 1 jp_end = cv_size_jb + cv_size_je + cv_size_jp jl_start = cv_size_jb + cv_size_je + cv_size_jp + 1 jl_end = cv_size_jb + cv_size_je + cv_size_jp + cv_size_jl jt_start = jl_end + 1 jt_end = jl_end + cv_size_jt grad_jo = 0.0 grad_jb = 0.0 grad_je = 0.0 grad_jd = 0.0 grad_jp = 0.0 grad_js = 0.0 grad_jl = 0.0 grad_jm = 0.0 grad_jt = 0.0 inc_div = 0.0 jcdf_flag = .false. !------------------------------------------------------------------------- ! [1.0] calculate grad_v (jo): !------------------------------------------------------------------------- call da_allocate_y(iv, jo_grad_y) #if (WRF_CHEM == 1) call da_allocate_y_chem_sfc(iv, jo_grad_y) #endif if (present(re)) then call da_calculate_grady(iv, re, jo_grad_y) if ( iter > 0 .and. test_gradient ) jcdf_flag = .true. #ifdef VAR4D call da_transform_vtoy_adj(cv_size, be, grid%ep, grad_jo, iv, & grid%vp, grid%vv, xbx, jo_grad_y, grid, config_flags, jcdf_flag, grid%vp6, grid%vv6) #else call da_transform_vtoy_adj(cv_size, be, grid%ep, grad_jo, iv, & grid%vp, grid%vv, xbx, jo_grad_y, grid, config_flags, jcdf_flag & !!When vchem includes initial conditions, need to add vchem=grid%vchem at end of this call #if (WRF_CHEM == 1) , vchem=grid%vchem & #endif ) #endif else #ifdef VAR4D call da_transform_vtoy(cv_size, be, grid%ep, cv, iv, grid%vp, & grid%vv, xbx, y, grid, config_flags, grid%vp6, grid%vv6) #else call da_transform_vtoy(cv_size, be, grid%ep, cv, iv, grid%vp, & grid%vv, xbx, y, grid, config_flags & #if (WRF_CHEM == 1) , vchem=grid%vchem & #endif ) #endif call da_calculate_grady(iv, y, jo_grad_y) #ifdef VAR4D call da_transform_vtoy_adj(cv_size, be, grid%ep, grad_jo, iv, & grid%vp, grid%vv, xbx, jo_grad_y, grid, config_flags, .true., grid%vp6, grid%vv6) #else call da_transform_vtoy_adj(cv_size, be, grid%ep, grad_jo, iv, & grid%vp, grid%vv, xbx, jo_grad_y, grid, config_flags, .true. & #if (WRF_CHEM == 1) , vchem=grid%vchem & #endif ) #endif grad_jo = - grad_jo !! Compensate for sign in calculation of grad_v (Jo) end if call da_deallocate_y(jo_grad_y) !------------------------------------------------------------------------- ! [2.0] calculate grad_v (jb): !------------------------------------------------------------------------- if (cv_size_jb > 0) grad_jb(1:cv_size_jb) = jb_factor * cv(1:cv_size_jb) !------------------------------------------------------------------------- ! [3.0] calculate grad_v (je): !------------------------------------------------------------------------- if (cv_size_je > 0) grad_je(je_start:je_end) = je_factor * cv(je_start:je_end) !---------------------------------------------------------------------- ! [3.1] calculate grad_v (jd): !---------------------------------------------------------------------- if (use_wpec) then if (var4d) call da_error(__FILE__,__LINE__,(/'Cannot use 4dvar with dynamic constraint'/)) if (wpec_factor <= 0) call da_error(__FILE__,__LINE__,(/'"wpec_factor" for dynamic constraint must be greater than zero'/)) grid%xa%grad_p_x(:,:,:)=0.0 grid%xa%grad_p_y(:,:,:)=0.0 call da_transform_vtod_wpec(cv_size, be, grid%ep, cv, grid%vp, grid%vv, xbx, grid) grid%xa%grad_p_x=(grid%xa%grad_p_x)/wpec_factor grid%xa%grad_p_y=(grid%xa%grad_p_y)/wpec_factor call da_transform_vtod_wpec_adj(cv_size, be, grid%ep, grad_jd, grid%vp, grid%vv, xbx, grid) end if !------------------------------------------------------------------------- ! [4.0] calculate grad_v (jp): !------------------------------------------------------------------------- #if defined(RTTOV) || defined(CRTM) if (use_varbc .and. cv_size_jp > 0) then do inst = 1, iv % num_inst do ichan = 1, iv%instid(inst)%nchan npred = iv%instid(inst)%varbc(ichan)%npred if (npred <= 0) cycle !! VarBC channels only do ipred = 1, npred id = iv%instid(inst)%varbc(ichan)%index(ipred) bgerr = iv%instid(inst)%varbc(ichan)%bgerr(ipred) if (bgerr > 0.0) & grad_jp(id) = (1/sqrt(bgerr)) * & SUM(cv(id) * iv%instid(inst)%varbc(ichan)%vtox(ipred,1:npred)) end do end do end do end if #endif !------------------------------------------------------------------------- ! [5.0] calculate grad_v (js): !------------------------------------------------------------------------- if (ANY(use_satcv)) then do inst = 1, iv % num_inst do n = iv%instid(inst)%info%n1, iv%instid(inst)%info%n2 ! loop for pixel ! Skin Temperature !----------------- if (use_satcv(1)) & grad_js(iv%instid(inst)%cv_index(n)%ts) = cv(iv%instid(inst)%cv_index(n)%ts) ! Cloud cover(s) !--------------- if (use_satcv(2)) then grad_js(iv%instid(inst)%cv_index(n)%cc) = cv(iv%instid(inst)%cv_index(n)%cc) WHERE (cv(iv%instid(inst)%cv_index(n)%cc) < 0.0 .or. & cv(iv%instid(inst)%cv_index(n)%cc) > 1.0 ) & grad_js(iv%instid(inst)%cv_index(n)%cc) = grad_js(iv%instid(inst)%cv_index(n)%cc) + & 10.0 * cv(iv%instid(inst)%cv_index(n)%cc) end if end do end do end if !------------------------------------------------------------------------- ! [6.0] calculate grad_v (jl): !------------------------------------------------------------------------- if (cv_size_jl > 0) grad_jl(jl_start:jl_end) = cv(jl_start:jl_end) !------------------------------------------------------------------------- ! [6.1] calculate grad_v (jm): !------------------------------------------------------------------------- if (use_divc) then call da_transform_vtox(grid, cv_size, xbx, be, grid%ep, cv, grid%vv, grid%vp) if ( be%ne > 0 .and. alphacv_method == alphacv_method_xa ) then call da_transform_vpatox(grid, be%ne, grid%ep, grid%vp) call da_add_xa(grid%xa, grid%xa_ens) !grid%xa = grid%xa + xa_ens end if call da_transform_xtoxa(grid) call da_divergence_constraint(grid%xb, grid%xa%u, grid%xa%v, inc_div) inc_div = inc_div/divc_factor call da_zero_x(grid%xa) call da_divergence_constraint_adj(grid, grid%xa%u, grid%xa%v, inc_div) call da_transform_xtoxa_adj(grid) if (be % ne > 0 .and. alphacv_method == alphacv_method_xa) then call da_transform_vpatox_adj(grid, be%ne, grid%ep, grid%vp) end if call da_zero_vp_type(grid%vp) call da_transform_vtox_adj(grid, cv_size, xbx, be, grid%ep, grid%vp, grid%vv, grad_jm) end if !------------------------------------------------------------------------- ! [n.n] calculate grad_t (jt): !------------------------------------------------------------------------- if (cv_size_jt > 0) then do n = 1, iv%varbc_tamdar%nair do iphase = 1, iv%varbc_tamdar%nphase if (iv%varbc_tamdar%nobs(iphase,n) > 0 .and. iv%varbc_tamdar%ifuse(iphase,n) > 0) then do ipred = 1, iv%varbc_tamdar%npred id = iv%varbc_tamdar%index(ipred,iphase,n) bgerr = iv%varbc_tamdar%bgerr(ipred,iphase,n) if (bgerr > 0.0) & grad_jt(id) = (1/sqrt(bgerr)) * SUM(cv(id) * & iv%varbc_tamdar%vtox(ipred,1:iv%varbc_tamdar%npred,iphase,n)) end do end if end do end do end if !-------------------------------------------------------------------------------------------------- ! [7.0] calculate grad_v (j) = grad_v (jb) + grad_v (jo) + grad_v (je) + grad_v (jd) + grad_v (jp) + grad_v (js) + grad_v (jl) !-------------------------------------------------------------------------------------------------- grad = grad_jo + grad_jb + grad_je + grad_jd + grad_jp + grad_js + grad_jl + grad_jm + grad_jt !------------------------------------------------------------------------- ! [8.0] write Gradient: !------------------------------------------------------------------------- if (rootproc) then if (it == 1 .and. iter == 0) then write(unit=grad_unit,fmt='(a)')'Outer EPS Inner G Gb Go Ge Gd Gp Gs Gl Gm Gt' write(unit=grad_unit,fmt='(a)')'Iter Iter ' end if end if if ( present(re) ) then #if (WRF_CHEM != 1) gnorm_j = sqrt(da_dot_cv(cv_size, grad, grad, grid, be%cv_mz, be%ncv_mz, jp_start, jp_end)) gnorm_jo = sqrt(da_dot_cv(cv_size, grad_jo, grad_jo, grid, be%cv_mz, be%ncv_mz)) gnorm_jb = sqrt(da_dot_cv(cv_size, grad_jb, grad_jb, grid, be%cv_mz, be%ncv_mz)) gnorm_je = sqrt(da_dot_cv(cv_size, grad_je, grad_je, grid, be%cv_mz, be%ncv_mz)) gnorm_jd = sqrt(da_dot_cv(cv_size, grad_jd, grad_jd, grid, be%cv_mz, be%ncv_mz)) gnorm_jp = sqrt(da_dot_cv(cv_size, grad_jp, grad_jp, grid, be%cv_mz, be%ncv_mz, jp_start, jp_end)) gnorm_js = sqrt(da_dot_cv(cv_size, grad_js, grad_js, grid, be%cv_mz, be%ncv_mz)) gnorm_jl = sqrt(da_dot_cv(cv_size, grad_jl, grad_jl, grid, be%cv_mz, be%ncv_mz)) gnorm_jm = sqrt(da_dot_cv(cv_size, grad_jm, grad_jm, grid, be%cv_mz, be%ncv_mz)) gnorm_jt = sqrt(da_dot_cv(cv_size, grad_jt, grad_jt, grid, be%cv_mz, be%ncv_mz)) #else gnorm_j = sqrt(da_dot_cv(cv_size, grad, grad, grid, be%cv_mz, be%ncv_mz,be%cv_mz_chemic, jp_start, jp_end)) gnorm_jo = sqrt(da_dot_cv(cv_size, grad_jo, grad_jo, grid, be%cv_mz, be%ncv_mz,be%cv_mz_chemic)) gnorm_jb = sqrt(da_dot_cv(cv_size, grad_jb, grad_jb, grid, be%cv_mz, be%ncv_mz,be%cv_mz_chemic)) gnorm_je = sqrt(da_dot_cv(cv_size, grad_je, grad_je, grid, be%cv_mz, be%ncv_mz,be%cv_mz_chemic)) gnorm_jd = sqrt(da_dot_cv(cv_size, grad_jd, grad_jd, grid, be%cv_mz, be%ncv_mz,be%cv_mz_chemic)) gnorm_jp = sqrt(da_dot_cv(cv_size, grad_jp, grad_jp, grid, be%cv_mz, be%ncv_mz,be%cv_mz_chemic, jp_start, jp_end)) gnorm_js = sqrt(da_dot_cv(cv_size, grad_js, grad_js, grid, be%cv_mz, be%ncv_mz,be%cv_mz_chemic)) gnorm_jl = sqrt(da_dot_cv(cv_size, grad_jl, grad_jl, grid, be%cv_mz, be%ncv_mz,be%cv_mz_chemic)) #endif if (rootproc) & write(grad_unit,fmt='(2x,i2,1x,e10.3,2x,i4,10(1x,f10.3))') & it, eps(it), iter, gnorm_j, gnorm_jb, gnorm_jo, gnorm_je, gnorm_jd, gnorm_jp, gnorm_js, gnorm_jl, gnorm_jm, gnorm_jt end if if (trace_use) call da_trace_exit("da_calculate_gradj") end subroutine da_calculate_gradj