LCOV - code coverage report
Current view: | top level - dynamics/se - apply_iop_forcing.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 62 | 0.0 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 0 | 2 | 0.0 %

              Line data    Source code
           1             : module apply_iop_forcing_mod
           2             :
           3             : use shr_kind_mod,   only:r8 => shr_kind_r8, i8 => shr_kind_i8
           4             : use pmgrid,         only:plev, plevp, plon
           5             : use constituents,   only:pcnst, cnst_get_ind, cnst_name
           6             : use physconst,      only:rair,cpair
           7             : use cam_logfile,    only:iulog
           8             : use hybvcoord_mod,  only: hvcoord_t
           9             : use scamMod,        only: use_3dfrc, single_column, have_u, have_v, divT3d, divq3d, divt, divq, &
          10             :                           wfld, uobs, vobs, tobs, qobs, plevs0, have_divt3d, have_divq3d, &
          11             :                           scm_relax_bot_p,scm_relax_linear,scm_relax_tau_bot_sec, &
          12             :                           scm_relax_tau_sec,scm_relax_tau_top_sec,scm_relax_top_p, &
          13             :                           scm_relaxation,scm_relax_fincl,qinitobs
          14             :
          15             : use cam_abortutils, only: endrun
          16             : use string_utils,   only: to_upper
          17             :
          18             : implicit none
          19             :
          20             : public advance_iop_forcing
          21             : public advance_iop_nudging
          22             :
          23             : !=========================================================================
          24             : contains
          25             : !=========================================================================
          26             :
          27           0 : subroutine advance_iop_forcing(scm_dt, ps_in, &                   ! In
          28             :                     u_in, v_in, t_in, q_in, t_phys_frc, q_phys_frc, hvcoord, & ! In
          29             :                     u_update, v_update, t_update, q_update)       ! Out
          30             :
          31             : !-----------------------------------------------------------------------
          32             : !
          33             : ! Purpose:
          34             : ! Apply large scale forcing for t, q, u, and v as provided by the
          35             : !   case IOP forcing file.
          36             : !
          37             : ! Author:
          38             : ! Original version: Adopted from CAM3.5/CAM5
          39             : ! Updated version for E3SM: Peter Bogenschutz (bogenschutz1@llnl.gov)
          40             : !  and replaces the forecast.F90 routine in CAM3.5/CAM5/CAM6/E3SMv1/E3SMv2
          41             : !
          42             : !-----------------------------------------------------------------------
          43             :
          44             :   ! Input arguments
          45             :   real(r8), intent(in) :: ps_in             ! surface pressure [Pa]
          46             :   real(r8), intent(in) :: u_in(plev)        ! zonal wind [m/s]
          47             :   real(r8), intent(in) :: v_in(plev)        ! meridional wind [m/s]
          48             :   real(r8), intent(in) :: t_in(plev)        ! temperature [K]
          49             :   real(r8), intent(in) :: q_in(plev,pcnst)  ! q tracer array [units vary] already vertically advected
          50             :   real(r8), intent(in) :: t_phys_frc(plev)  ! temperature forcing from physics [K/s]
          51             :   real(r8), intent(in) :: q_phys_frc(plev,pcnst)  ! change in q due to physics.
          52             :   type (hvcoord_t), intent(in)   :: hvcoord
          53             :   real(r8), intent(in) :: scm_dt            ! model time step [s]
          54             :
          55             :   ! Output arguments
          56             :   real(r8), intent(out) :: t_update(plev)      ! updated temperature [K]
          57             :   real(r8), intent(out) :: q_update(plev,pcnst)! updated q tracer array [units vary]
          58             :   real(r8), intent(out) :: u_update(plev)      ! updated zonal wind [m/s]
          59             :   real(r8), intent(out) :: v_update(plev)      ! updated meridional wind [m/s]
          60             :
          61             :   ! Local variables
          62             :   real(r8) pmidm1(plev)  ! pressure at model levels
          63             :   real(r8) pintm1(plevp) ! pressure at model interfaces
          64             :   real(r8) pdelm1(plev)  ! pdel(k)   = pint  (k+1)-pint  (k)
          65             :   real(r8) t_lsf(plev)       ! storage for temperature large scale forcing
          66             :   real(r8) q_lsf(plev,pcnst) ! storage for moisture large scale forcing
          67             :   real(r8) fac, t_expan
          68             :
          69             :   integer i,k,m           ! longitude, level, constituent indices
          70             :
          71             :   character(len=*), parameter :: subname = 'advance_iop_forcing'
          72             :
          73             :   ! Get vertical level profiles
          74           0 :   call plevs0(plev, ps_in, hvcoord%ps0, hvcoord%hyam, hvcoord%hybm, hvcoord%hyai, hvcoord%hybi, pintm1 ,pmidm1 ,pdelm1)
          75             :
          76             :   !  Advance T and Q due to large scale forcing
          77           0 :   if (use_3dfrc) then
          78           0 :     if(.not.(have_divt3d.and.have_divq3d)) call endrun(subname//": FATAL: divt3d and divq3d not available")
          79           0 :     t_lsf(:plev) = divt3d(:plev)
          80           0 :     q_lsf(:plev,:pcnst) = divq3d(:plev,:pcnst)
          81             :   else
          82           0 :     t_lsf(:plev) = divt(:plev)
          83           0 :     q_lsf(:plev,:pcnst) = divq(:plev,:pcnst)
          84             :   endif
          85             :
          86           0 :   do k=1,plev
          87             :      ! Initialize thermal expansion term to zero.  This term is only
          88             :      !  considered if three dimensional forcing is not provided by IOP forcing file.
          89           0 :      t_expan = 0._r8
          90             :
          91           0 :     if (.not. use_3dfrc) then
          92           0 :       t_expan = scm_dt*wfld(k)*t_in(k)*rair/(cpair*pmidm1(k))
          93             :     endif
          94             :
          95           0 :      if (use_3dfrc) then
          96           0 :         do m=1,pcnst
          97             :            ! When using 3d dynamics tendencies, SCM skips the vertical advection step and thus
          98             :            ! q_in at this point has not had physics tendencies applied
          99           0 :            q_update(k,m) = q_in(k,m) + scm_dt*(q_phys_frc(k,m) + q_lsf(k,m))
         100             :         end do
         101           0 :         t_update(k)   = t_in(k) + t_expan + scm_dt*(t_phys_frc(k) + t_lsf(k))
         102             :      else
         103           0 :         do m=1,pcnst
         104             :            ! When not using 3d dynamics tendencies, q_in at this point has had physics tend
         105             :            ! applied and has been vertically advected. Only horizontal dyn tend needed for forecast.
         106           0 :            q_update(k,m) = q_in(k,m) + scm_dt*q_lsf(k,m)
         107             :         end do
         108           0 :         t_update(k)   = t_in(k) + t_expan + scm_dt*t_lsf(k)
         109             :      end if
         110             :   end do
         111             :
         112             :   !  Set U and V fields
         113             :
         114           0 :   if ( have_v .and. have_u ) then
         115           0 :     do k=1,plev
         116           0 :       u_update(k) = uobs(k)
         117           0 :       v_update(k) = vobs(k)
         118             :     enddo
         119             :   endif
         120             :
         121           0 : end subroutine advance_iop_forcing
         122             :
         123             : !=========================================================================
         124             :
         125           0 : subroutine advance_iop_nudging(ztodt, ps_in,                        &      ! In
         126             :                                tfcst, qfcst, ufcst, vfcst, hvcoord, &      ! Inout
         127             :                                relaxt, relaxq )                            ! Out
         128             :
         129             :   !-----------------------------------------------------------------------
         130             :   !
         131             :   ! Purpose:
         132             :   ! Option to nudge t and q to observations as specified by the IOP file
         133             :   !-----------------------------------------------------------------------
         134             :
         135             :   ! Input arguments
         136             :   real(r8), intent(in) :: ztodt            ! model time step [s]
         137             :   real(r8), intent(in) :: ps_in            ! surface pressure [Pa]
         138             :   type (hvcoord_t), intent(in)   :: hvcoord
         139             :
         140             :   ! Output arguments
         141             :   real(r8), intent(inout) :: tfcst(plev)  ! updated temperature [K]
         142             :   real(r8), intent(inout) :: qfcst(plon,plev,pcnst)  ! updated const field
         143             :   real(r8), intent(inout) :: ufcst(plev)  ! updated U wind
         144             :   real(r8), intent(inout) :: vfcst(plev)  ! updated V wind
         145             :   real(r8), intent(out) :: relaxt(plev)    ! relaxation of temperature [K/s]
         146             :   real(r8), intent(out) :: relaxq(plev)    ! relaxation of vapor [kg/kg/s]
         147             :
         148             :   ! Local variables
         149             :   integer :: i, k, m
         150             :   real(r8) pmidm1(plev)  ! pressure at model levels
         151             :   real(r8) pintm1(plevp) ! pressure at model interfaces
         152             :   real(r8) pdelm1(plev)  ! pdel(k)   = pint  (k+1)-pint  (k)
         153             :
         154             :   ! --------------------------- !
         155             :   ! For 'scm_relaxation' switch !
         156             :   ! --------------------------- !
         157             :
         158             :   real(r8) rtau(plev)
         159             :   real(r8) relax_T(plev)
         160             :   real(r8) relax_u(plev)
         161             :   real(r8) relax_v(plev)
         162             :   real(r8) relax_q(plev,pcnst)
         163             :   ! +++BPM: allow linear relaxation profile
         164             :   real(r8) rslope ! [optional] slope for linear relaxation profile
         165             :   real(r8) rycept ! [optional] y-intercept for linear relaxtion profile
         166             :   logical scm_fincl_empty
         167             :
         168             :   ! ------------------------------------------------------------------- !
         169             :   ! Relaxation to the observed or specified state                       !
         170             :   ! We should specify relaxation time scale ( rtau ) and                !
         171             :   ! target-relaxation state ( in the current case, either 'obs' or 0 )  !
         172             :   ! ------------------------------------------------------------------- !
         173             :
         174           0 :   if ( .not. scm_relaxation) return
         175             :
         176           0 :   call plevs0(plev, ps_in, hvcoord%ps0, hvcoord%hyam, hvcoord%hybm, hvcoord%hyai, hvcoord%hybi, pintm1 ,pmidm1 ,pdelm1)
         177             :
         178           0 :   relax_T(:)             = 0._r8
         179           0 :   relax_u(:)             = 0._r8
         180           0 :   relax_v(:)             = 0._r8
         181           0 :   relax_q(:plev,:pcnst)  = 0._r8
         182             :   ! +++BPM: allow linear relaxation profile
         183             :   ! scm_relaxation is a logical from scamMod
         184             :   ! scm_relax_tau_top_sec and scm_relax_tau_bot_sec are the relaxation times at top and bottom of layer
         185             :   ! also defined in scamMod
         186           0 :   if ( scm_relax_linear ) then
         187           0 :      rslope = (scm_relax_top_p - scm_relax_bot_p)/(scm_relax_tau_top_sec - scm_relax_tau_bot_sec)
         188           0 :      rycept = scm_relax_tau_top_sec - (rslope*scm_relax_top_p)
         189             :   endif
         190             :
         191           0 :   scm_fincl_empty=.true.
         192           0 :   do i=1,pcnst
         193           0 :      if (len_trim(scm_relax_fincl(i)) > 0) then
         194           0 :         scm_fincl_empty=.false.
         195           0 :         scm_relax_fincl(i)=trim(to_upper(scm_relax_fincl(i)))
         196             :      end if
         197             :   end do
         198             :
         199           0 :   do k = 1, plev
         200           0 :      if ( pmidm1(k) <= scm_relax_bot_p.and.pmidm1(k) >= scm_relax_top_p ) then ! inside layer
         201           0 :         if (scm_relax_linear) then
         202           0 :            rtau(k) = rslope*pmidm1(k) + rycept ! linear regime
         203             :         else
         204           0 :            rtau(k)         = max( ztodt, scm_relax_tau_sec ) ! constant for whole layer / no relax outside
         205             :         endif
         206           0 :      else if  (scm_relax_linear .and. pmidm1(k) <= scm_relax_top_p ) then ! not linear => do nothing / linear => use upper value
         207           0 :         rtau(k) = scm_relax_tau_top_sec ! above layer keep rtau equal to the top
         208             :      endif
         209             :      ! +BPM: this can't be the best way...
         210             :      ! I put this in because if rtau doesn't get set above, then I don't want to do any relaxation in that layer.
         211             :      ! maybe the logic of this whole loop needs to be re-thinked.
         212           0 :      if (rtau(k) /= 0) then
         213           0 :         relax_T(k)      = -  ( tfcst(k)     - tobs(k) )    / rtau(k)
         214           0 :         relax_u(k)      = -  ( ufcst(k)     - uobs(k) )    / rtau(k)
         215           0 :         relax_v(k)      = -  ( vfcst(k)     - vobs(k) )    / rtau(k)
         216           0 :         relax_q(k,1)    = -  ( qfcst(1,k,1) - qobs(k) )    / rtau(k)
         217           0 :         do m = 2, pcnst
         218           0 :            relax_q(k,m) = -  ( qfcst(1,k,m) - qinitobs(k,m)   )    / rtau(k)
         219             :         enddo
         220           0 :         if (scm_fincl_empty .or. ANY(scm_relax_fincl(:) == 'T')) &
         221           0 :              tfcst(k)        =      tfcst(k)     + relax_T(k)   * ztodt
         222           0 :         if (scm_fincl_empty .or.ANY(scm_relax_fincl(:) == 'U')) &
         223           0 :              ufcst(k)        =      ufcst(k)     + relax_u(k)   * ztodt
         224           0 :         if (scm_fincl_empty .or. ANY(scm_relax_fincl(:) == 'V')) &
         225           0 :              vfcst(k)        =      vfcst(k)     + relax_v(k)   * ztodt
         226           0 :         do m = 1, pcnst
         227           0 :            if (scm_fincl_empty .or. ANY(scm_relax_fincl(:) == trim(to_upper(cnst_name(m)))) ) then
         228           0 :               qfcst(1,k,m) =      qfcst(1,k,m) + relax_q(k,m) * ztodt
         229             :            end if
         230             :         enddo
         231             :      end if
         232             :   enddo
         233             :
         234             : end subroutine advance_iop_nudging
         235             :
         236             : !-----------------------------------------------------------------------
         237             :
         238             : end module apply_iop_forcing_mod

Generated by: LCOV version 1.14